Thermal Analysis of a Functionally Graded Coating/Substrate System Using the Approximated Transfer Approach

As a heterogeneous material, functionally graded material (FGM) behaves as continuously changed material properties in certain directions from one composition to another, and hence it has received much more attention for biomedical applications and thermal protections to achieve innovative functions that conventional homogeneous material cannot accomplish. However, due to the particularly small thickness ratio of coating to substrate in practice, the conventional mesh discretization of the coating region is inefficient. To simplify the meshing procedure and increase the efficiency of analysis, the approximated transfer algorithm based on the concept of finite difference is developed for transferring boundary conditions applied on the coating surface to the interface of coating and substrate. As a result, only the substrate with transferred convection boundary conditions needs to be solved numerically, i.e., by the fundamental-solution based hybrid finite element method (HFS-FEM) with high accuracy and feasible polygonal element construction, in which only integrals along the element boundary are evaluated because of the application of fundamental solutions of the problem as kernel functions of interior approximated fields. Finally, numerical experiments including the single-layered, multi-layered and functionally graded coatings are carried out to verify the accuracy and applicability of the present method.


Introduction
Practically, biomaterials are usually not homogeneous and show certain functional gradation formed by biological adaptation [1].For instance, the porous bone is graded from the external cortical bone, a dense and stiff structure, to the internal cancellous bone which is a porous one [2].A similar feature can be observed in bamboo [3] and mollusc shell [4].Such a functionally graded structure optimizes the material's response to external loading.Therefore, it is of great importance that the artificial structures for implant have similar natural functional gradations.As an example of bone implants for knee joint replacement, a functionally graded interlayer is usually designed to improve the acceptance of artificial implants by living tissues.This can be achieved for titanium implants functionalized with a graded coating [5].Besides, in the case of dental implants, the functionally graded layers were introduced in the titanium and its alloys to increase osteoconductivity and biomechanical bonding between the implant and the surrounding bony tissue [2].As an application of the HA-TiO 2 -Ti functionally graded biocompatible coating system designed by Kumar and Wang for implant applications [6], the hydroxyapatite (HA) powder was mixed with the titanium oxide (TiO 2 ) power in different weight percentages to form five thin coating layers and then attached to the Ti6A14V metal substrate to obtain the resulting functionally graded composite.Besides, the multi-layered or functionally graded coating film was also used as a thermal barrier coating for protecting the covered substrate under severe thermal-mechanical environments [7][8][9][10].For such multi-layered or functionally graded coating/substrate systems, the prediction of temperature distribution in the covered substrate domain is critical for further analysis of thermal stresses and deformations to control them within acceptable levels.
To analyze the multi-layered coating system or the functionally graded coating (FGC) system, numerical modelling is widely employed, including finite element modelling [11,12] and boundary element modelling [13,14].However, it is observed that the thickness of coating is usually very thin, i.e., within µm scale, compared to the size of the covered substrate domain in millimeter scale.The primary obstacle to the development of efficient and accurate computational modelling for such a thin coating system is that the meshing operation in the thin coating layers becomes very difficult and accuracy and efficiency becomes worse as the thickness of the coating layer decreases dramatically.
To deal with such an obstacle, the transfer approximation based on finite difference technique [13] and the nonlinear transformation technique [14,15] were respectively developed for transferring the specified boundary conditions on the coating surface to the interface of the coating layer and substrate so that only the substrate region needs to be solved, instead of the whole coating/substrate system.However, there are a very few studies that focused on thermal analyses of heterogeneous coating systems, i.e., FGC [6,16,17].
This paper pays attention to the first attempt of the thermal analysis of a functionally graded coating/substrate system by combining the approximated transfer algorithm and the hybrid finite element with fundamental solution kernels, an effective numerical tool called fundamental-solution based hybrid finite element method (HFS-FEM), which can be viewed as the combination of finite element method (FEM) [12,18] and boundary element method (BEM) [19][20][21] and has such inherent advantages as the stiffness equation containing element boundary integrals only, specially-purposed elements, arbitrary polygonal-shaped element construction, and better accuracy over the conventional FEM.However, it is worth noting that the dependence of fundamental solutions of HFS-FEM may limit its application to those problems without the explicit expression of fundamental solutions.Detailed summaries of the development of HFS-FEM and the state-of-the-art can be found in [22,23].In this study, the approximated transfer algorithm based on finite difference is developed to establish a recurrence formula to treat the multi-layered coating and the functionally graded coating.The applied thermal boundary conditions on the top surface of coating can be transferred to the surface of substrate in convection form and then only the substrate domain is solved by the HFS-FEM to analyze the thermal behavior in it.Finally, the accuracy and applicability of the present method is discussed by several numerical examples involving the single-layered coating, the multi-layered coating and the functionally graded coating.

Problem Statement
To address the generalization of the present method, a general two-dimensional functionally graded coating/substrate system depicted in Figure 1 is considered in the study, where the thickness of heterogeneous coating and substrate is h c and h s , respectively.It is assumed that the substrate region is homogeneous and isotropic, while the coating layer is graded along the thickness of it.Besides, it is worth pointing out that the substrate is not limited to the finite domain shown in Figure 1, and the half-infinite substrate domain is also workable with the numerical procedure developed below.For the composite system consisting of the single-layered functionally graded coating and the substrate, the basic partial differential equations at point x = (x1,x2) governing the steady-state temperature distribution in the coating region Ωc and the substrate region Ωs can be expressed as: where qc and qs are heat flux vectors in the coating and substrate regions, respectively, and can be expressed in terms of temperature variable by the Fourier's law, In Equation (2), kc and ks are the thermal conductivity of the coating material and the substrate, respectively.Tc and Ts are the corresponding temperature fields in each region, respectively.For the homogeneous coating, the thermal conductivity kc is constant, while it may change in terms of spatial coordinates along the thickness of the functionally graded coating.In this study, the following two graded forms are considered, as shown in Figure 1, for the linear form ( ) for the exponental form bx a bx k x ae where a and b are constant coefficients, which can be determined by the continuous change of thermal conductivity from one material phase to another material phase.Substituting Equation (2) into Equation (1) yields, Obviously, the two governing equations defined in Equation ( 1) are independent of each other, owing to the material difference in thermal conductivity.To couple them, the temperature and normal heat flux continuous conditions at the interface Гi of the coating and substrate should be added for representing local thermal equilibrium: on the interface 0 on the interface where nc = −ns are respectively the unit outward normal vector to the interface for the coating domain and substrate domain.
Typically, the thin functionally graded material (FGM) coating is practically laminated by several homogeneous sublayers with different proportions along its thickness [6,16].For such cases, the governing Equations ( 4) can be rewritten in general form: For the composite system consisting of the single-layered functionally graded coating and the substrate, the basic partial differential equations at point x = (x 1 ,x 2 ) governing the steady-state temperature distribution in the coating region Ω c and the substrate region Ω s can be expressed as: where q c and q s are heat flux vectors in the coating and substrate regions, respectively, and can be expressed in terms of temperature variable by the Fourier's law, In Equation (2), k c and k s are the thermal conductivity of the coating material and the substrate, respectively.T c and T s are the corresponding temperature fields in each region, respectively.For the homogeneous coating, the thermal conductivity k c is constant, while it may change in terms of spatial coordinates along the thickness of the functionally graded coating.In this study, the following two graded forms are considered, as shown in Figure 1, a + bx 2 for the linear form ae bx 2 for the exponental form where a and b are constant coefficients, which can be determined by the continuous change of thermal conductivity from one material phase to another material phase.Substituting Equation (2) into Equation (1) yields, Obviously, the two governing equations defined in Equation ( 1) are independent of each other, owing to the material difference in thermal conductivity.To couple them, the temperature and normal heat flux continuous conditions at the interface Γ i of the coating and substrate should be added for representing local thermal equilibrium: where n c = −n s are respectively the unit outward normal vector to the interface for the coating domain and substrate domain.
Typically, the thin functionally graded material (FGM) coating is practically laminated by several homogeneous sublayers with different proportions along its thickness [6,16].For such cases, the governing Equations (4) can be rewritten in general form: where k cm is thermal conductivity of the m th layer, T cm denotes the temperature in the m th layer and M is the number of homogeneous sublayers in the functionally graded coating domain, as shown in Figure 2. Note that in Equation ( 6), the (M + 1) th layer refers to the substrate domain and Subsequently, the continuous conditions at the interface of adjacent layers m and m + 1 can be written as: T cm = T c(m+1) q cm n cm + q c(m+1) n c(m+1) = 0 (7) where q cm = −k cm ∆T cm is the heat flux vector in the m th (m = 1, 2, • • • , M + 1) layer, and n cm = −n c(m+1) is the normal vector to the interface.
Coatings 2018, 8, x FOR PEER REVIEW 4 of 17 where kcm is thermal conductivity of the m th layer, Tcm denotes the temperature in the m th layer and M is the number of homogeneous sublayers in the functionally graded coating domain, as shown in Figure 2. Note that in Equation ( 6), the (M + 1) th layer refers to the substrate domain and kc(M+1) = ks, Tc(M+1) = Ts, Ωc(M+1) = Ωs.
Subsequently, the continuous conditions at the interface of adjacent layers m and m + 1 can be written as:

Approximated Transfer Algorithm for Thin Graded Coating Layer
Conventional multi-region numerical methods such as the FEM and the HFS-FEM are inefficient when the coating is particularly thin.Under such ultra-thin circumstances, extremely refined mesh with massive integration points is required.To tackle this problem, the approximated transfer approach for the single homogeneous coating [13] is firstly reviewed and then is extended to deal with the multi-layered coating and the functionally graded coating by an iterative scheme.Although the temperature distribution in the thin coating domain cannot be given through the present approximated transfer algorithm, the algorithm can significantly improve the computational efficiency in the substrate region, which is particularly interesting to assess the "bioactive" and "protective" effect of the coating.

Thin Single-Layered Homogeneous Coating
For the single-layer homogeneous coating shown in Figure 3, it is assumed that a mixed boundary condition is applied at point Qc, ( ) where α and β0 are constant coefficients, and γ is the given value, respectively.Practically, the combinations α = 1 and β0 = 0, α = 0 and β0 = 1, and α ≠ 0 and β0 ≠ 0 respectively correspond to the temperature condition, heat flux condition, and convection condition.
According to the forward finite difference technique [24], the first-derivative of the temperature field Tc at point Qs can be approximated as:

Approximated Transfer Algorithm for Thin Graded Coating Layer
Conventional multi-region numerical methods such as the FEM and the HFS-FEM are inefficient when the coating is particularly thin.Under such ultra-thin circumstances, extremely refined mesh with massive integration points is required.To tackle this problem, the approximated transfer approach for the single homogeneous coating [13] is firstly reviewed and then is extended to deal with the multi-layered coating and the functionally graded coating by an iterative scheme.Although the temperature distribution in the thin coating domain cannot be given through the present approximated transfer algorithm, the algorithm can significantly improve the computational efficiency in the substrate region, which is particularly interesting to assess the "bioactive" and "protective" effect of the coating.

Thin Single-Layered Homogeneous Coating
For the single-layer homogeneous coating shown in Figure 3, it is assumed that a mixed boundary condition is applied at point Q c , where α and β 0 are constant coefficients, and γ is the given value, respectively.Practically, the combinations α = 1 and β 0 = 0, α = 0 and β 0 = 1, and α = 0 and β 0 = 0 respectively correspond to the temperature condition, heat flux condition, and convection condition.
According to the forward finite difference technique [24], the first-derivative of the temperature field T c at point Q s can be approximated as: from which the normal heat flux at point Q s can be given by: Introducing the interfacial conditions (5) we have: from which such relation can be derived as: Besides, the heat flux term in Equation ( 8) is assumed to be same as that applied at the coating substrate interface [13], that is, Then, substituting Equations ( 12) and ( 13) into the mixed boundary condition ( 8) produces: from which it is found that the applied boundary condition (8) on the top surface of the coating is approximately transferred to the interface of substrate and coating.Obviously, the results of the specified temperature condition which corresponds to α = 1, β 0 = 0 and the specified heat flux condition which corresponds to α = 0, β 0 = 1 can be regarded as special cases of the general form (14).
Coatings 2018, 8, x FOR PEER REVIEW 5 of 17 from which the normal heat flux at point Qs can be given by: Introducing the interfacial conditions (5) we have: from which such relation can be derived as: Besides, the heat flux term in Equation ( 8) is assumed to be same as that applied at the coating substrate interface [13], that is, Then, substituting Equations ( 12) and ( 13) into the mixed boundary condition (8) produces: from which it is found that the applied boundary condition (8) on the top surface of the coating is approximately transferred to the interface of substrate and coating.Obviously, the results of the specified temperature condition which corresponds to α = 1, β0 = 0 and the specified heat flux condition which corresponds to α = 0, β0 = 1 can be regarded as special cases of the general form (14).

Thin Functionally Graded Coating
Once the approximated transfer approach for the single-layered homogeneous coating is derived, the procedure can be extended for deriving the recurrence formula for the functionally graded coating.In this study, the functionally graded coating is generally divided into M thin layers, and each layer is approximately homogeneous, as shown in Figure 2. It is assumed that the thickness and thermal conductivity of the m th layer is hcm and kcm If a general form of boundary condition is applied at the top surface point Qc1 of the layer #1,

Thin Functionally Graded Coating
Once the approximated transfer approach for the single-layered homogeneous coating is derived, the procedure can be extended for deriving the recurrence formula for the functionally graded coating.In this study, the functionally graded coating is generally divided into M thin layers, and each layer is approximately homogeneous, as shown in Figure 2. It is assumed that the thickness and thermal conductivity of the m th layer is h cm and k cm (m = 1, 2, • • • , M), respectively.
by means of the derivation of the single-layered case, the transferred boundary condition at point Q c2 on the interface of layer #1 and layer #2 can be written as: with Following this transferring rule, the resulting thermal-response at the point Q s on the upper surface of the substrate domain can be written as: in which the transfer coefficient β M is given by the following recurrence relation: Specially, when the temperature condition is applied on the coating surface, that is, α = 1, β 0 = 0, Equation ( 18) reduces to If the normal heat flux condition is given on the coating surface, that is, α = 0, β 0 = 1, Equation (18) reduces to (q s n s ) Q s = γ .Besides, the procedure above shows that the present approximated transfer approach is also suitable for multi-layered homogeneous or inhomogeneous coating.

Hybrid Element Formulation
From the above transfer procedure, it is found that the substrate region should be solved under the approximated convection-typed boundary condition (18) applied on the coating-substrate interface.This can be done through the HFS-FEM including the convection term, which is implemented by introducing fundamental solution kernels for element interior fields and shape functions for element boundary frame field.Compared to the conventional FEM and BEM, and the Trefftz FEM [25,26], this method has some inherent characteristics like stiffness equations containing element boundary integrals only, specially-purposed elements and polygonal-shaped elements and has been successfully applied to elastic [22,27,28] and thermal conduction [23,[29][30][31] problems with better accuracy over the conventional FEM and more flexible multi-sided element construction, prior to this study.Additionally, different to the BEM, the HFS-FEM can be applied to multi-domain or multi-material problems more conveniently.
Before deriving the computing formulation, the mixed boundary condition given in Equation ( 18) should be reorganized to meet the requirement for implementing the HFS-FEM.In the hybrid variational functional of HFS-FEM for thermal analysis, three types of boundary conditions are involved: specified temperature condition, specified heat flux condition and specified convection condition.For such cases, Equation ( 18) is rewritten as the following convection form: where respectively denote the generalized convection coefficient and environment temperature.For two-dimensional steady-state heat transfer in the isotropic homogeneous substrate media, the weak form of the double-variable hybrid functional for a particular hybrid finite element e, i.e., element #1 in Figure 4, can be expressed by [29].(22) in which q represents the specified value of normal heat flux q n = q s •n on the element heat flux boundary Γ qe , n = (n 1 ,n 2 ) is the unit normal vector to the element boundary, T s and T s are respectively the temperature within the element and on the element boundary, Ω e and Γ e respectively denote the element domain and element boundary, and Γ ce and Γ qe stand for the element convection boundary and heat flux boundary, respectively.It is noted that the four integral terms on the right side of Equation ( 22) represent the effects of thermal energy equilibrium, specified normal heat flux condition, inter-element continuity condition, and the convection condition, respectively.Besides, it is assumed that the coordinates x 1 and x 2 are set along the element plane.The temperature in the element is only a function of two space coordinates x 1 and x 2 , when the prescribed temperature and heat flux are independent of the thickness direction of element.This means that the kinematics of heat flow along the thickness direction of the element can be ignored.Therefore, the two-dimensional hybrid finite element with thickness t(x 1 , x 2 ) can be used to model the computing domain.Here, the constant element thickness is considered in this work, for simplicity.For two-dimensional steady-state heat transfer in the isotropic homogeneous substrate media, the weak form of the double-variable hybrid functional for a particular hybrid finite element e, i.e., element #1 in Figure 4, can be expressed by [29].In the application of variational functional (22), the temperature and its gradient fields at any point x in the interior of the element  are approximated by linear combination of fundamental solutions centered at series of source points s s ( 1, , )  In the application of variational functional (22), the temperature and its gradient fields at any point x in the interior of the element e are approximated by linear combination of fundamental solutions centered at series of source points x s j (j = 1, • • • , n s ) locating on the pseudo boundary similar to the element boundary Γ e , that is, where c e = {c ej } T (j = 1, • • • , n s ) is the coefficient vector consisting of unknown source intensities at n s source points locating outside the element domain.N and T are respectively coefficient matrices consisting of temperature fundamental solutions T*(x,x j s ) and heat flux fundamental solutions Coatings 2019, 9, 51 8 of 17 For the two-dimensional isotropic and homogeneous substrate taken into consideration, the temperature fundamental solutions are given as [19].
Because the interior fields are independently defined for adjacent elements, i.e., element #1 and #2 in Figure 4, they should be weakly linked through the common interface.In order to enforce the conformity on the common interface of adjacent elements, the element frame temperature field T s can be defined in terms of the generalized nodal temperature vector where N denotes the matrix consisting of the conventional one-dimensional shape functions widely used in the standard FEM [18] and BEM [19], n d is the number of nodes of the element e.
The substitution of the interior and frame fields related to the element e into the functional (22) gives where and It is worth pointing out that during the procedure of deriving Equation (29) involving integrals along the element boundary only, the distinct feature that the interior fields exactly satisfy the governing equation is used to remove the domain integral in the functional (22).Moreover, the element boundary integrals in (29) mean that the polygonal element with more sides can be flexibly constructed for practical computation, not limited to triangle and quadrilateral elements in the conventional FEM.
Finally, the stationary conditions of Π me with respect to the coefficient vector c e and the nodal temperature vector d e , yield the optional relationship between c e and d e , and the element stiffness equation, where It is clear from the above procedure that the element stiffness matrix is symmetric, and after assembling the global stiffness matrix keeps sparse and symmetric features, such as that in the conventional FEM.Additionally, the evaluation of the right-handed vector (33) is the same as that in the conventional FEM.Therefore, these features make the implementation of the present hybrid finite element conveniently incorporate into the existing FEM program.

Numerical Results
In this section, five tests including the single-layered homogeneous coating with simple and complicated boundary conditions and the functionally graded coatings with exponential and linear graded forms are carried out for demonstrating the accuracy and applicability of the present method.Since the focus of the current work is on the development of generalized computational technique for analyzing various coating/substrate systems, the outmost material phase of the FGC layer or the single-layered homogeneous coating is arbitrarily assumed to have a thermal conductivity 6 W m −1 K −1 .The thermal conductivity of the substrate material is assumed to have a value of 28 W m −1 K −1 .These values can be adjusted to meet the practical coating/substrate system.Additionally, the convergence of the present hybrid finite element and the effect of the location of source points have been well investigated for thermal analysis (see literature [23,[29][30][31] for reference), so they are not discussed in this study.Moreover, the 4-sided hybrid finite element with 4 edges and 8 nodes is employed without exception in the following computation.

Homogeneous Thin Coating with Simple Boundary Conditions
Firstly, a simple test with a simple temperature boundary condition on the coating surface is performed to investigate the accuracy of the present method.In this test shown in Figure 5, it is assumed that both the left and right sides of the coating/substrate system are insulated, and the top and bottom sides are of specified temperature f c and f s .From the basic heat transfer theory, the one-dimensional analytical solutions of temperature for the coating and substrate domains are given by from which it is seen that temperature distributions both in the coating and in the substrate are linear.
In the practical computation, it is assumed that the coating/substrate model is composed of a square substrate with h s = 1 mm and a thin homogeneous coating with thickness h c .The specified temperature boundary condition on the top and bottom surfaces of the model is assumed to be f c = 1173 K and f s = 298 K, respectively.For the purpose of comparing the numerical and analytical results, the temperature distributions are evaluated in the substrate domain for different coating thickness values, which are controlled by the specific thickness ratio h c /h s of the coating and substrate ranging from 10 −6 to 10 −1 , and in each analysis, the hybrid finite element mesh with an element size of 0.1 mm is identically modelled.Since only the substrate domain is solved, the present approximated model is clearly more computationally efficient than the full-domain model consisting of the coating and substrate domains.
For different coating thickness values, the determined temperature values at three specific points respectively locating on the interface and within the substrate region, i.e., A (0.4 mm, 1 mm), B (0.2 mm, 0.9 mm) and C (0.6 mm, 0.4 mm), are recorded for analysis.Results in Table 1 indicate that there is perfect agreement between the numerical results and the analytical solutions, even for very small coating thicknesses.Thus, the present method shows a significant potential for analyzing thin coating problems.Moreover, it is seen that the temperature is nearly independent of the thickness ratio h c /h s when its value is equal to or less than 10 −4 .It is reasonable that the extremely thin coating will make the temperature profile in the coating layer very close to the specified value f c .
However, it is worth noting that the test above is a one-dimensional linear case and it is not enough to assess the capability of the present method.Hence, more complex two-dimensional problems should be considered to check whether the present method can be used for them.coating problems.Moreover, it is seen that the temperature is nearly independent of the thickness ratio hc/hs when its value is equal to or less than 10 −4 .It is reasonable that the extremely thin coating will make the temperature profile in the coating layer very close to the specified value fc.
However, it is worth noting that the test above is a one-dimensional linear case and it is not enough to assess the capability of the present method.Hence, more complex two-dimensional problems should be considered to check whether the present method can be used for them.

Homogeneous Thin Coating with Complicated Boundary Conditions
In the second test, a complicated two-dimensional heat transfer is tested.The same geometrical domain and element mesh as those in the first test are taken here.Along the outer boundary of the coating/substrate system in Figure 5, the temperature boundary conditions are applied by the following exact solutions of temperature respectively for the coating and substrate regions where η is any constant.It is obvious that Equation ( 36) is constructed by the polynomial including 2 nd and 1 st terms to satisfy the governing Equation ( 1) and the interfacial condition ( 5) so that the two-dimensional heat transfer in the coating/substrate system is ensured.

Homogeneous Thin Coating with Complicated Boundary Conditions
In the second test, a complicated two-dimensional heat transfer is tested.The same geometrical domain and element mesh as those in the first test are taken here.Along the outer boundary of the coating/substrate system in Figure 5, the temperature boundary conditions are applied by the following exact solutions of temperature respectively for the coating and substrate regions where η is any constant.It is obvious that Equation ( 36) is constructed by the polynomial including 2nd and 1st terms to satisfy the governing Equation ( 1) and the interfacial condition (5) so that the two-dimensional heat transfer in the coating/substrate system is ensured.Firstly, we take η = 10 in Equation ( 36) such that the analytical temperature distribution is independent of the thickness ratio h c /h s .Results at three different points the same as those in the first test are listed in Table 2, from which it is found that the accuracy of numerical results begins to slightly decrease as the thickness ratio drops below 10 −2 , and the relative error to the analytical solution at the interfacial point A is the largest (0.32%), compared to that at the interior points B and C, i.e., 0.21% at the point B and 0.05% at the point C. It is reasonable that the points B and C are far away from the interface on which the approximated transfer boundary conditions are applied.Besides, Figure 6 shows the variation of temperature along the interface and we observe that, although the temperature changes nonlinearly in terms of the coordinate x 1 , the numerical results calculated by the present method are in good agreement with the exact solutions.
Subsequently, we take η = 1000 h c /h s so that the analytical temperature fields both in the coating and in the substitute severely depend on the thickness ratio h c /h s .Results in Table 3 indicate that the numerical results from the present method are still in good agreement with the exact solutions, although the values of temperature at the three specific points increase with the increase of the thickness Coatings 2019, 9, 51 ratio.Besides, the temperature variations along the interface for different thickness ratios are plotted in Figure 7 and the good agreement between the numerical results and the analytical solutions is observed again.although the temperature changes nonlinearly in terms of the coordinate x1, the numerical results calculated by the present method are in good agreement with the exact solutions.Subsequently, we take η = 1000 hc/hs so that the analytical temperature fields both in the coating and in the substitute severely depend on the thickness ratio hc/hs.Results in Table 3 indicate that the numerical results from the present method are still in good agreement with the exact solutions, although the values of temperature at the three specific points increase with the increase of the thickness ratio.Besides, the temperature variations along the interface for different thickness ratios are plotted in Figure 7 and the good agreement between the numerical results and the analytical solutions is observed again.Table 3. Temperature results at the three specific points for different thickness ratios for η = 1000 hc/hs in the second test.

Partially Heated Homogeneous Thin Coating
A substrate with a partially heated homogeneous thin coating is considered in the third test.The boundary conditions specified in Figure 8a are complicated enough to ensure two-dimensional heat transfer.The side length ℎ of the square substrate domain is 7 mm, and the local length  of

Partially Heated Homogeneous Thin Coating
A substrate with a partially heated homogeneous thin coating is considered in the third test.The boundary conditions specified in Figure 8a are complicated enough to ensure two-dimensional heat transfer.The side length h s of the square substrate domain is 7 mm, and the local length l 0 of the specific extremely high temperature constraint is 1.75 mm.The related hybrid finite element mesh in the substrate domain is generated with 700 elements, as shown in Figure 8b.

Partially Heated Homogeneous Thin Coating
A substrate with a partially heated homogeneous thin coating is considered in the third test.The boundary conditions specified in Figure 8a are complicated enough to ensure two-dimensional heat transfer.The side length ℎ of the square substrate domain is 7 mm, and the local length  of the specific extremely high temperature constraint is 1.75 mm.The related hybrid finite element mesh in the substrate domain is generated with 700 elements, as shown in Figure 8b.Because there is no analytical solution for this test, the numerical results from the present method are compared to the available BEM solutions [13].To do so, the determined temperature values at the point (1 mm, 6.41 mm) within the substrate region are recorded in Table 4 for comparison under different coating thicknesses between 1 µm to 1 mm.From Table 4, it is shown that the deviation between numerical results from the single-domain BEM and the present method can be ignored.Thus, the present method can be used for engineering analysis of thin coating.Besides, it is seen from Table 4 that, when the coating thickness is less than 10 µm, the temperature at the specific point just changes slightly.This can be attributed to the extremely small thickness of the coating layer.Because there is no analytical solution for this test, the numerical results from the present method are compared to the available BEM solutions [13].To do so, the determined temperature values at the point (1 mm, 6.41 mm) within the substrate region are recorded in Table 4 for comparison under different coating thicknesses between 1 µm to 1 mm.From Table 4, it is shown that the deviation between numerical results from the single-domain BEM and the present method can be ignored.Thus, the present method can be used for engineering analysis of thin coating.Besides, it is seen from Table 4 that, when the coating thickness is less than 10 µm, the temperature at the specific point just changes slightly.This can be attributed to the extremely small thickness of the coating layer.

Linear-Form Functionally Graded Coating
Different to the conventional homogeneous coating studied above, the FGC offers the possibility of significantly improving interface reliability caused by material mismatch [17], and thus is considered in this and the next tests.Identical geometry and boundary conditions to that depicted in Figure 5 are employed.It is assumed that the thermal conductivity of the FGC varies along the thickness of coating in the following linear form: in which the coefficients a and b can be determined from the continuous change of thermal conductivity from the outmost material phase k c0 = 6 W m −1 K −1 to the inner substrate material phase k S = 28 W m −1 K −1 , i.e., when x 2 = h c + h s , k c = k c0 , and when x 2 = h s , k c = k s , from which we have From Equation (38), it is seen that the value of a is linearly dependent of the thickness ratio h s /h c , but b depends on the coating thickness only.For example, for h c = 0.001 mm, we have a = 22,028 and b = −22,000.
From the basic heat transfer theory, the one-dimensional analytical solution of temperature in the coating/substrate system can be derived by substituting Equation (37) into Equation (4), To demonstrate the capability of the present method for analyzing this inhomogeneous coating, the coating is uniformly divided into 2, 6 and 10 thin layers, respectively, and each layer is assumed to be homogeneous.The related thermal conductivity of each layer is evaluated by the coordinate of centroid of the layer.Table 5 demonstrates the temperature results at the interfacial point (0.5 mm, 1 mm) for different coating thicknesses.It is clearly seen that the numerical accuracy of the present method becomes better as the number of layers increases for each value of the coating thickness.Besides, for M = 2, the results deteriorate quickly as the coating thickness is less than 0.1 mm; however, more layers, i.e., M = 6 and M = 10, can produce very consistent results with the exact solutions even though the thickness is larger than 0.1 mm.Finally, the temperature variations along the middle line x 1 = 0.5 mm for different coating thicknesses are plotted in Figure 9 when the 10-layered coating model is employed.It is found that the temperature profile in the substrate tends to be unchanged as the thickness of the coating becomes small.

Exponential-Form Functionally Graded Coating
In the last test, the same model as that in the fourth test is considered and the exponent form, From Equation (41), it is seen that the value of a is seriously dependent on the thickness ratio hc/hs, while the coefficient b depends on the coating thickness only.Too small hc/hs will cause extremely large a, which is not suitable for practical use.For example, for the case of hc = 0.1 mm, hs = 1 mm, we have a = 1.3716 × 10 8 , b = −15.4045.
From the basic heat transfer theory, the one-dimensional analytical solution of temperature in the coating and substrate domains can be derived by substituting Equation (40) into Equation (4)

, b h h b h h b h h bh ax b h h b h h bh b h h b h h bh e k f f abe f h f k e f k e T e abe h k e e abe h k e e
To demonstrate the performance of the present method for exponentially graded coating, the coating is also divided into 2, 6, 10 thin layers and each layer is assumed to be homogeneous.When the coating thickness changes from 0.1 to 0.3 mm, the temperature results calculated from the

Exponential-Form Functionally Graded Coating
In the last test, the same model as that in the fourth test is considered and the exponent form, is utilized to represent the continuous change of the thermal conductivity of the coating material.
The coefficients a and b are similarly determined from the continuous change of thermal conductivity from the outmost material phase k c0 = 6 W m −1 K −1 to the inner substrate material phase k s = 28 W m −1 K −1 , i.e., when x 2 = h c + h s , k c = k c0 , and when x 2 = h s , k c = k s , from which we have From Equation (41), it is seen that the value of a is seriously dependent on the thickness ratio h c /h s , while the coefficient b depends on the coating thickness only.Too small h c /h s will cause extremely large a, which is not suitable for practical use.For example, for the case of h c = 0.1 mm, h s = 1 mm, we have a = 1.3716 × 10 8 , b = −15.4045.
From the basic heat transfer theory, the one-dimensional analytical solution of temperature in the coating and substrate domains can be derived by substituting Equation (40) into Equation ( 4) To demonstrate the performance of the present method for exponentially graded coating, the coating is also divided into 2, 6, 10 thin layers and each layer is assumed to be homogeneous.When the coating thickness changes from 0.1 to 0.3 mm, the temperature results calculated from the present method both in Table 6 and in Figure 10 show good agreement with the analytical solutions.Moreover, the convergent tendency of temperature can be observed from Table 6 when the number of layers increases, as expected.Additionally, in contrast to the results in Table 5, it is seen that the exponentially graded distribution can produce a smaller temperature at the interface than the linear distribution under the same coating thickness.This is because the exponential variation of thermal conductivity can yield a smaller value than the linear form.the number of layers increases, as expected.Additionally, in contrast to the results in Table 5, it is seen that the exponentially graded distribution can produce a smaller temperature at the interface than the linear distribution under the same coating thickness.This is because the exponential variation of thermal conductivity can yield a smaller value than the linear form.

Conclusions
In this work, the general functionally graded coating/substrate system is solved by integrating the approximated transfer algorithm and the HFS-FEM for predicting the temperature distribution in the substrate and numerical results, demonstrating that the proposed method can be applied to various coatings with high computational efficiency and accuracy.In contrast to the previous works in the literature, the proposed method can be summarized as follows:

•
The approximated recurrence formula based on the finite difference technique is derived for approximately transferring the general thermal boundary condition applied on the coating surface to the coating/substrate interface.

•
The established transfer approach is highly suitable for treating thin multi-layered coating or functionally graded coating.

•
With the present transfer approach, only the substrate needs to be solved, instead of the whole coating/substrate system.

•
With the fundamental solutions of the problem as kernel functions, the HFS-FEM involving the evaluation of element boundary integrals only is implemented in the substrate region to achieve high accuracy and feasible element construction.

Conclusions
In this work, the general functionally graded coating/substrate system is solved by integrating the approximated transfer algorithm and the HFS-FEM for predicting the temperature distribution in the substrate and numerical results, demonstrating that the proposed method can be applied to various coatings with high computational efficiency and accuracy.In contrast to the previous works in the literature, the proposed method can be summarized as follows: • The approximated recurrence formula based on the finite difference technique is derived for approximately transferring the general thermal boundary condition applied on the coating surface to the coating/substrate interface.

•
The established transfer approach is highly suitable for treating thin multi-layered coating or functionally graded coating.

•
With the present transfer approach, only the substrate needs to be solved, instead of the whole coating/substrate system.

•
With the fundamental solutions of the problem as kernel functions, the HFS-FEM involving the evaluation of element boundary integrals only is implemented in the substrate region to achieve high accuracy and feasible element construction.

Figure 1 .
Figure 1.Configuration of the functionally graded coating/substrate system.

Figure 1 .
Figure 1.Configuration of the functionally graded coating/substrate system.

Figure 2 .
Figure 2. Schematic diagrams of (a) the multi-layered coating/substrate system, and (b) the adjacent layers.

Figure 2 .
Figure 2. Schematic diagrams of (a) the multi-layered coating/substrate system, and (b) the adjacent layers.

Figure 3 .
Figure 3. Configuration of the single-layered homogeneous coating/substrate system.

Figure 3 .
Figure 3. Configuration of the single-layered homogeneous coating/substrate system.
form of boundary condition is applied at the top surface point Q c1 of the layer #1, 22) in which  represents the specified value of normal heat flux qn = qs•n on the element heat flux boundary Гqe, n = (n1,n2) is the unit normal vector to the element boundary, Ts and  are respectively the temperature within the element and on the element boundary, Ωe and Гe respectively denote the element domain and element boundary, and Гce and Гqe stand for the element convection boundary and heat flux boundary, respectively.It is noted that the four integral terms on the right side of Equation (22) represent the effects of thermal energy equilibrium, specified normal heat flux condition, inter-element continuity condition, and the convection condition, respectively.Besides, it is assumed that the coordinates  and  are set along the element plane.The temperature in the element is only a function of two space coordinates  and  , when the prescribed temperature and heat flux are independent of the thickness direction of element.This means that the kinematics of heat flow along the thickness direction of the element can be ignored.Therefore, the twodimensional hybrid finite element with thickness ( ,  ) can be used to model the computing domain.Here, the constant element thickness is considered in this work, for simplicity.

Figure 4 .
Figure 4. Schematic diagram of two-dimensional hybrid finite element.
on the pseudo boundary similar to the element boundary Гe, that is, ce = {cej} T (j = 1, ••• , ns) is the coefficient vector consisting of unknown source intensities at ns source points locating outside the element domain.N and T are respectively coefficient matrices consisting of temperature fundamental solutions T*(x,xj s ) and heat flux fundamental solutions

Figure 4 .
Figure 4. Schematic diagram of two-dimensional hybrid finite element.

Figure 5 .
Figure 5. Configuration of the coating/substrate system in the first test.

Figure 5 .
Figure 5. Configuration of the coating/substrate system in the first test.

Figure 6 .
Figure 6.Variation of temperature along the interface for η = 10 in the second test.

Figure 6 .
Figure 6.Variation of temperature along the interface for η = 10 in the second test.

Figure 7 .
Figure 7. Variation of temperature along the interface for η = 1000 hc/hs in the second test.

Figure 7 .
Figure 7. Variation of temperature along the interface for η = 1000 h c /h s in the second test.

Figure 7 .
Figure 7. Variation of temperature along the interface for η = 1000 hc/hs in the second test.

Figure 8 .
Figure 8.(a) Configuration of the coating/substrate system in the third test, (b) the in-plane hybrid finite element mesh in the substrate domain.

Figure 8 .
Figure 8.(a) Configuration of the coating/substrate system in the third test, (b) the in-plane hybrid finite element mesh in the substrate domain.

Figure 9 .
Figure 9. Temperature variation along the middle line  = 0.5 mm for different coating thicknesses for linear gradation.
is utilized to represent the continuous change of the thermal conductivity of the coating material.The coefficients a and b are similarly determined from the continuous change of thermal conductivity from the outmost material phase kc0 = 6 W m −1 K −1 to the inner substrate material phase ks = 28 W m −1 K −1 , i.e., when x2 = hc + hs, kc = kc0, and when x2 = hs, kc = ks, from which we have

Figure 9 .
Figure 9. Temperature variation along the middle line x 1 = 0.5mm for different coating thicknesses for linear gradation.

Figure 10 .
Figure 10.Temperature variation along the middle line  = 0.5 mm for different coating thicknesses (10 layers) for exponential gradation

Funding:
The research was funded by the National Natural Science Foundation of China (No. 11472099), the fund of the State Key Laboratory of Structural Analysis for Industrial Equipment in the Dalian University of Technology (No. GZ1610) and the Program for Innovative Research Team of Science & Technology of Henan Province of China (No. 19IRTSTHN020).

Figure 10 .
Figure 10.Temperature variation along the middle line x 1 = 0.5 mm for different coating thicknesses (10 layers) for exponential gradation.
e c e − d e T g e + c e T G e d e −

Table 2 .
Temperature results at the three specific points for different thickness ratios for η = 10 in the second test.

Table 2 .
Temperature results at the three specific points for different thickness ratios for η = 10 in the second test.

Table 3 .
Temperature results at the three specific points for different thickness ratios for η = 1000 h c /h s in the second test.h

c /h s Point A Point B Point C EXACT Present EXACT Present EXACT Present
Coatings 2018, 8, x FOR PEER REVIEW 12 of 17

Table 4 .
Results at the specific point for different coating thickness in the third test.

Table 4 .
Results at the specific point for different coating thickness in the third test.

Table 5 .
Temperature at the interfacial point (0.5 mm, 1 mm) for different coating thicknesses for linear gradation.

Table 6 .
Temperature at the interfacial point (0.5 mm, 1 mm) for different coating thicknesses for exponential gradation.

Table 6 .
Temperature at the interfacial point (0.5 mm, 1 mm) for different coating thicknesses for exponential gradation.