Surface contact behavior of functionally graded thermoelectric materials indented by a conducting punch

The contact problem for thermoelectric materials with functionally graded properties is considered. The material properties, such as the electric conductivity, the thermal conductivity, the shear modulus, and the thermal expansion coefficient, vary in an exponential function. Using the Fourier transform technique, the electro-thermoelastic problems are transformed into three sets of singular integral equations which are solved numerically in terms of the unknown normal electric current density, the normal energy flux, and the contact pressure. Meanwhile, the complex homogeneous solutions of the displacement fields caused by the gradient parameters are simplified with the help of Euler’s formula. After addressing the non-linearity excited by thermoelectric effects, the particular solutions of the displacement fields can be assessed. The effects of various combinations of material gradient parameters and thermoelectric loads on the contact behaviors of thermoelectric materials are presented. The results give a deep insight into the contact damage mechanism of functionally graded thermoelectric materials (FGTEMs).


Introduction
The temperature difference between cold and hot extremities of thermoelectric devices can be directly converted into electrical energy by the Seebeck effect. Conversely, through the Peltier effect, the thermoelectric device can be energized to achieve cooling [1][2] . Therefore, they have broad applications in the fields of aviation detectors, engine waste heat power generation, thermoelectric refrigeration, and wearable equipment [3][4] .
conducting and electrically conducting conductor and has the total electric current J e0 , the total energy flux J u0 , and the external force P without any loss being permitted within the contact area. Fig. 1 The geometry of the plane contact problem for the FGTEMs (color online)

Basic equations
The thermo-electric-elastic properties of the FGTEMs are approximated by γ(y) = γ 0 e δy , λ(y) = λ 0 e χy , G(y) = G 0 e αy , β(y) = β 0 e μy (1) for 0 < y < ∞, where γ(y), λ(y), G(y), and β(y) are the electric conductivity, the thermal conductivity, the shear modulus, and the thermal expansion coefficient, respectively. The corresponding values at the surface y = 0 of the FGTEMs are denoted as γ 0 , λ 0 , G 0 , and β 0 . The notations δ, χ, α, and μ are the corresponding gradient parameters. These gradient parameters can be either positive or negative. When the gradient index is positive, the value of the corresponding material property approaches infinity as the depth increases, while for negative gradient indexes, it approaches zero. The arbitrary material properties of FGTEMs may be more realistic. As for classic functionally graded materials, the arbitrary parameters have been investigated [22][23]25] . The present exponential model of FGTEMs may serve as the solid foundation for the future work of FGTEMs with a further arbitrary variation of properties.
With the parameters introduced above, by considering that both electric charges and energy are conserved, the governing equations and the constitutive equations can be written as ∇ · j e = 0, ∇ · j u = 0, (2) j e = −γ(y)∇V − γ(y)ε∇T, q = εT j e − λ(y)∇T, where ∇· = ∂ ∂x + ∂ ∂y is a divergence operator. ∇ = ∂ ∂x , ∂ ∂y T is a gradient operator. j e = (j ex , j ey ) T denotes the electric current density. q = (q x , q y ) T represents the thermal flux. The superscript "T" means transpose of a vector or a matrix. The energy flux j u can be expressed as j u = q + j e V . ε, V , and T are the Seebeck coefficient, the electric potential, and the temperature, respectively. To facilitate the solution procedure, the electrochemical potential function F is defined as F = V + εT . After calculating, the governing equations and the constitutive equations in Eqs. (2) and (3) revealing the thermoelectric gradient properties can be rewritten as follows: The electrochemical potential function F (x, y) can be solved from the first equation of Eq. (4). The temperature T (x, y) can be obtained from the second equation of Eq. (4) which can further help to solve the elastic field.
The thermoelasticity constitutive relations of the FGTEMs in terms of the temperature T (x, y) as well as the displacement components u(x, y) and w(x, y) are written as where κ = 3 − 4ν, and β * = (1 + v)β 0 for the considered plane strain case. The equilibrium equations of the FGTEMs are given as follows: To solve the governing equations in Eqs. (4) and (7), the classical Fourier integral transform technique and the superposition principle are adopted.

Electric field
The solution to the first equation of Eq. (4) under the regularity condition at infinity (F (x, ∞) = 0, |x| < ∞) can be obtained as follows: where n(s) is negative and given as n(s) = − δ 2 − δ 2 4 + s 2 , and A(s) is to be determined. Due to the problem considered here, the surface normal electric current density j ey (x, 0) can be denoted as where J * e (x) is the unknown normal electric current density under the punch, and satisfies in which J e0 represents the total electric current acted on the punch. Considering the first equation of Eq. (5) and Eqs. (8) and (9), A(s) can be given as (8) and performing an asymptotic analysis, one can get where the generalized Fredholm integral kernel K 11 (x, r) without singularity is Furthermore, by considering the boundary condition F (x, 0) = F 0 , |x| < a, Eq.
where J * e (r) satisfies the equilibrium equation (10). It is noted that K 11 (x, r) will vanish when the electric conductivity gradient parameter δ equals zero. Then, Eq. (13) becomes one that describes the isotropic thermoelectric materials with the following exact solution [29] : J * e (x) = Je0 π √ a 2 −x 2 , |x| < a, which can verify our derivation.

Temperature field
The distributions of the temperature function T (x, y) and the normal energy flux j uy (x, y) on the surface or at infinity can be expressed as where T 0 and J * u (x) are the constant value temperature within the contact area and the unknown normal energy flux inside the contact region, respectively. The prescribed total energy flux J u0 above the punch is Indeed, the gradient index of electric conductivity δ and that of thermal conductivity χ for real FGTEMs may be different, which will make it hard to obtain the particular solution of the temperature function in an analytical form when solving the second equation of Eq. (4). Therefore, the assumption that δ and χ have the same values is used to obtain the particular solution of temperature fields in an analytical form. Except for other statements, it is modeled that χ is equal to δ in the following calculation. Now, by employing the method used in the electric field and considering the second equation of Eq. (4) and boundary conditions in Eqs. (14) and (15), the temperature field involving the particular solutions and homogeneous solutions can be obtained as follows: where B(s) is given as with the equilibrium equation (16) to complete the solution.
The above assumption δ = χ is used to obtain the particular solution of temperature fields in an analytical form, which finally leads to the similar equations, i.e., Eqs. (13) and (18). While the total electric current J e0 and total energy flux J u0 appearing in the equilibrium equations (10) and (16) may make the two similar singular integral equations have different solutions, the results obtained here are primary and can act as the beginning of a more complex situation.

Mechanical field
Similarly, the Fourier integral transform technique can be used to solve Eq. (7). By considering the displacement boundary conditions at infinity, u(x, ∞) = 0, w(x, ∞) = 0, |x| < ∞, the complex homogeneous solutions to Eq. (7) in the Fourier transform field can be written as where , in which m j (j = 1, 2) with negative real parts are given by It is noticed that the two characteristic roots in Eq. (20) are conjugate complex roots caused by the FGTEMs. Then, m j (j = 1, 2) can be denoted as [30] where the bar represents conjugation. o 1 = Re(m 1 ), and τ 1 = Im(m 1 ). Re(·) and Im(·) are the real part and the imaginary part, respectively. By substituting Eq. (21) into Eq. (19) and using Euler's formula, the transformed complex homogeneous solutions can be rewritten in the real homogeneous solutions' form as follows: where Γ 1 = Re(M 1 ), and Θ 1 = Im(M 1 ).
With the convolution formula, the nonlinear term in Eq. (17) can be addressed. Then, the particular solutions to Eq. (7) under the thermoelectric effects are where η 1 (s, ξ) = n(ξ) + n(s − ξ), η 2 (s) = n(s), and E j , Z j , and Δ j (j = 1, 2) are given by , it can be seen that the parts relating to functions A and B represent the effects of the electric field and the temperature field on the displacement field, respectively.
As mentioned before, the gradient parameters of the material properties can be either positive or negative. These two different cases make the kernels of the integral equations sufficiently different. Whether the gradient parameters are positive or negative, the root n(s) is almost always negative, and the characteristic roots m 1 = m 2 given in Eq. (20) have negative real parts verified by the numerical tests. Thus, for the half-plane problem with y > 0, the displacements and the stresses vanish at infinity by checking Eqs. (19), (22), (23), and (6).

Solutions of the contact stress
With the transform matrix method, the relationship between the surface displacement and the surface stress can be obtained. Interested readers can find detailed information in Ref. [29]. The normal displacement can then be given as follows: where In the above equations, Differentiating Eq. (25) with respect to x yields The asymptotic behavior of the integrands for the variable s can be given as where . The surface normal stress inside the contact region is unknown, denoted as p(x), while it remains free outside the contact region, i.e., σ yy (x, 0) = −p(x), |x| < a; σ yy (x, 0) = 0, |x| > a.
The mechanical boundary conditions applied to Eq. (7) can be written as By separating the leading term in Eq. (28) and considering Eqs. (30) and (31), one can get the singular integral equation corresponding to the contact pressure as follows: where the generalized Fredholm integral kernel function K 22 (x, r) without singularity is and the right-hand side function f (x) revealing the thermoelectric effects and the gradient characteristic of thermoelectric materials is given as The obtained right-hand side expression of Eq. (32), which is first given, reveals the noteworthy difference between the contact analysis of classical materials (e.g., piezoelectric materials [31] and piezomagnetic materials [32] ) and that of FGTEMs considered here. Specifically, from Eq. (32) and Eqs. (34)-(36), one can see that the equivalent force f (x) consists of thermoelectric loads relating to the terms of A(·) in the electric field, B(·) in the temperature field, and the gradient parameters α and μ.
To obtain a physically possible solution, the following equilibrium condition should be satisfied: As for homogeneous thermoelectric materials, there exists sQ 21 = Q ∞ 1 which can lead to K 22 (x, r) = 0. Then, the singular integral equation (32) for FGTEMs can degenerate into the case of isotropic thermoelectric materials. Furthermore, it should be noted that if without thermoelectric loads and gradient parameters, L 2 (s) associated with the particular solutions of the displacements will also vanish. Then, the singular integral equation for isotropic materials without thermoelectric effects can be reduced to with the aid of Eq. (37), and the exact solution can be obtained as p(x) = P π √ a 2 −x 2 , |x| < a, which can illustrate the correctness of our derivation again.

Intensity factors
The unknown normal electric current density J * e (r), the unknown normal energy flux J * u (r), and the unknown contact pressure p(r) all satisfy the first kind Cauchy-type singular integral equation demonstrated in Eqs. (13), (18), and (32), which can be solved numerically with the collocation method. In the normalized interval as x = a x and r = a r, the functions J * e (r), J * u (r), and p(r) are where T j ( r) = cos(j arccos( r)) is the Chebyshev polynomial of the first kind.
With the help of the Gauss-Chebyshev quadrature method, one can obtain systems of linear algebraic equations as follows: where r n = cos 2n−1 2N π , n = 1, 2, · · · , N, and x l = cos l N π , l = 1, 2, · · · , N − 1. Once the coefficients Y j , U j , and B j are determined, the various intensity factors, such as the electric current density intensity factors K e (±a), the energy flux intensity factors K u (±a), and the stress intensity factors K I (±a), can be evaluated as follows: When the values of the above-calculated intensity factors reach a certain threshold, contact damage will occur due to crack initiation and propagation.

Numerical results and discussion
From the solution process, one can see that if without the gradient parameters δ, χ, α, and μ, the FGTEMs considered here will become homogeneous ones.
The numerical analysis is carried out to investigate the singularities near the contact edges including K e (±a), K u (±a), and K I (±a). The distributions of the contact pressure, the normal energy flux, and the normal electric current density below the rigid flat punch are also studied. The material parameters of the Bi 2 Te 3 -based thermoelectric material in Refs. [33] and [34] are selected as the surface properties required for this paper. Specifically, the electrical conductivity, the thermal conductivity, the thermal expansion coefficient, and the shear modulus at the surface y = 0 are taken as γ 0 = 1 × 10 5 S/m, λ 0 = 2.2 W/(m·K), β 0 = 1.68 × 10 −5 K −1 , and G 0 = 16.786 GPa, respectively.

Convergence studies
The convergence analyses of K 11 in Eq. (12) and K 22 in Eq. (33) are given in Figs. 2(a) and 2(b), respectively. The variation of the integral value of K 11 (a x l , a r n ) with the truncating number N t1 under different gradient indexes aδ = aχ = 1/4, 1/2, 2, 4 is presented in Fig. 2(a). One can see that the magnitudes of K 11 become steady as N t1 increases under aδ = aχ = 1/4, 1/2, while the magnitudes of K 11 have a slight fluctuation with aδ and aχ equaling 2 or 4. Therefore, the specific values of K 11 and their relative errors are listed in Table 1, which shows that the relative errors maintain within 4% when N t1 reaches 1 × 10 8 and are less than 2% once N t1 is beyond 2 × 10 8 . Thus, the value of N t1 is taken as 2 × 10 8 in the following calculations.  The kernel function K 22 in Eq. (33) is relevant to the gradient parameter of the shear modulus aα. Figure 2(b) demonstrates the effects of the magnitudes of the truncating number N t2 on the convergent behavior of K 22 when aα = 0.001, 1/3, 1, 3. When aα equals 3, the curve fluctuates slightly, and the corresponding relative errors are lower than 2% as N t2 is over 1.5×10 8 . Taking the remaining results of Fig. 2(b) into account, one may take N t2 as 2 × 10 8 in the following calculations, which can ensure the accuracy for the selected gradient parameters.
It should be pointed out that the efficiency of the collocation method sufficiently depends on the value of the gradient index. Table 2 gives the effects of the collocation number N with aδ = aχ = 1/4, 1/2, 2, 4 on the maximum normal electric current density j ey (0, 0)/j ey0 and the maximum normal energy flux j uy (0, 0)/j uy0 with j ey0 = J e0 /(2a) and j uy0 = J u0 /(2a). It can be seen that the results are closely equal to each other with larger N and the data with N = 35 or N = 50 are nearly identical. Thereby, N = 35 is employed in the thermoelectric calculations. Table 3 shows the effects of the electric conductivity gradient parameter aδ = 1/4, 1, 4, the thermal conductivity gradient parameter aχ = 1/4, 1, 4, the shear modulus gradient parameter aα = 0.001, 1, 2, and the thermal expansion coefficient gradient parameter aμ = 1, 2, 4 on the convergence behavior of the collocation method in the calculation of the maximum contact pressure σ yy (0, 0)/σ 0 with σ 0 = P/(2a). The effects of the collocation number on the maximum contact pressure σ yy (0, 0)/σ 0 under different values of the thermoelectric loads J u0 and J e0 are presented in Table 4. Consider all data given in Tables 3 and 4, and the results become stable with the increase in the collocation number under various gradient parameters. N = 25 or N = 30 can ensure sufficient accuracy. Thus, N = 30 is used in the following calculations of the contact pressure. Consider the four specific gradient parameters taken in Table 3, and the ranges of the gradient parameters for which the accuracy is acceptable are selected as aδ, aχ, aα, aμ 4.   Figure 3 reveals the effects of various gradient parameters aδ and aχ fixed at 1/4, 1/2, 0, 0.001, 2, and 4 on the normalized electric current density j ey (x, 0)/j ey0 and the normalized energy flux j uy (x, 0)/j uy0 . The closed-form results, which can be obtained as

Effects of the gradient parameters on thermoelectric fields
are validated with small values of the gradient parameters aδ = aχ= 0.001 as the degenerated results shown in Fig. 3. Excellent agreement can be found between the reduced results of aδ = aχ = 0.001 computed by the program and those of aδ = aχ = 0 obtained by the closed-form solutions, which further validates the present derivations. Furthermore, there are singularities in j ey (x, 0)/j ey0 and j uy (x, 0)/j uy0 at both ends of the contact region for the flat punch as shown in Fig. 3. It can be seen that when aδ and aχ take different magnitudes, the ranges of j ey (x, 0)/j ey0 and j uy (x, 0)/j uy0 generated at the contact surface also change. Therefore, by considering the actual operating environment of thermoelectric materials, the appropriate gradient parameters can be selected according to the theoretical analysis of this paper for the rational preparation of thermoelectric materials.  Table 5 shows the effects of different values of the gradient parameters aδ and aχ on the intensity factors K e (±a)/K 0 and K u (±a)/K 0 . K u (±a)/K 0 is equal to K e (±a)/K 0 due to the assumptions of aχ = aδ and J e0 = J u0 . Table 5 demonstrates that K e (±a)/K 0 and K u (±a)/K 0 decrease gradually as the gradient parameters increase. This also proves that the rational design of the gradient parameters of the thermoelectric materials can effectively reduce the singularity of the normal electric current density and the normal energy flux near the contact edges.

Effects of the gradient parameters and thermoelectric loads on the elastic field
The effects of aδ and aχ on σ yy (x, 0)/σ 0 are shown in Fig. 4 with aα= 1, aμ = 0, J e0 = 0.002 A/m, and J u0 = 0.002 W/m. Figure 4 shows that the gradient parameters aδ and aχ have no significant effects on the contact pressure. This shows that aδ and aχ of thermoelectric materials have limited effects on improving the elastic contact behavior of these materials.
As pointed out before, if without thermoelectric loads, i.e., J e0 = 0 A/m and J u0 = 0 W/m, the right-hand side function f (x) of the singular integral equation (32) is expected to vanish. This means that the distribution of contact pressure is only related to the shear modulus gradient parameter aα. When aα = 0, the closed-form solution of contact pressure can be obtained as p(x) = P/(π a 2 − x 2 ), |x| < a. Figure 5 shows the effects of aα on σ yy (x, 0)/σ 0 verified with the exact solution. It can be found that the degenerated results of aα = 0.001 evaluated by the program agree with those of aα = 0 given by the closed-form solutions, which again validates the present derivations. The effects of aα on σ yy (x, 0)/σ 0 are shown in Fig. 6, which depicts that as aα increases, σ yy (x, 0)/σ 0 in the middle part of the contact area significantly increases, while it decreases around the contact edges. Figure 6 also shows that when aα is small, the overall distribution of contact pressure is relatively uniform.  Figure 7 exhibits how aμ affects σ yy (x, 0)/σ 0 . It is clear that σ yy (x, 0)/σ 0 in the middle of the contact area slightly decreases as aμ increases. Figure 8 shows the effects of the total energy flux J u0 on σ yy (x, 0)/σ 0 . From Fig. 8, one can see that the increment of J u0 leads to the increase of σ yy (x, 0)/σ 0 in the middle of the contact area. Figure 9 depicts the effects of the total electric current J e0 on σ yy (x, 0)/σ 0 , and shows that as J e0 increases, the contact pressure in the middle of the contact area may be less compressive.

Conclusions
This paper proposes an analytical model of the FGTEMs indented by a rigid flat punch. The exponential variation law is used for the properties of the electric conductivity, the thermal conductivity, the shear modulus, and the thermal expansion coefficient. The convergent analysis of kernel functions and collocation methods under different gradient parameters is presented. Numerical tests are obtained to analyze the effects of these gradient parameters and the thermoelectric loads on the near-edge behavior of the FGTEMs. The following conclusions can be drawn.
(i) The electric conductivity gradient parameter and the thermal conductivity gradient parameter have significant effects on the distributions of the normal electric current density and the normal energy flux, while they have a slight effect on the distribution of the contact pressure.
(ii) Both a larger shear modulus gradient parameter and an enlarged total energy flux could lead to a higher contact pressure in the middle of the contact region.
(iii) As the thermal expansion coefficient gradient parameter and the total electric current increase, the contact pressure at the contact edge slightly increases.
(iv) The quantitative analysis in this paper should be helpful for the design of new FGTEMs and the further understanding of the contact behavior of FGTEMs.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.