Thermoelectric (TE) technology has attracted great attention because it can convert heat into electric power directly with the advantages of no moving parts and long durability, especially for the applications in deep space exploration, waste heat recovery, and other special fields^[1,2,3]. Over the past several decades, many high-performance TE materials and devices were developed^[4,5,6,7,8], providing great opportunities for large-scale application of TE power generation. Among them, skutterudite-based (SKD) devices exhibited very high conversion efficiencies up to 12%^[9], making SKD one of the most potential candidates for the practical applications.

As the SKD device steps forward to industrial application from laboratory^[10], the mechanical reliability becomes the top concern since TE devices often work under harsh service conditions, such as high temperature difference, mechanical vibration with wide frequency, extremely long service time. A TE device comprises of n-, p-type TE legs and electrodes as major components. Previously, researchers paid a lot of attention to the mechanical property of SKD materials themselves^[11,12,13], and demonstrated that the SKD materials exhibit excellent mechanical performance to bear service stress^[14,15]. However, at the electrode interface, the Sb elements of skutterudites matrix diffused into electrode materials such as Cu, Ni, and Mo^[16,17] or diffusion barrier layer (DBL) such as Ti, Mo, Mo-Ti, Cr₈₀Si₂₀, Nb, Ti_100-xAl_x (x= 3-12)^{[18,19,20,21,22,23]} during the long-time service at high temperatures. The elemental diffusions and chemical reactions at the interface not only result in dramatical increase of the interfacial electrical and thermal resistivity^[24], but also induce the mechanical injury or disability.

Numerical analysis for residual interfacial stress in multilayer system was well studied in the past decades. It was recognized that the coefficient of thermal expansion (CTE) of the component layers plays a critical role on interfacial reliability^[25]. Interface morphology should also affect interface stability in multilayer systems like thermal barrier coatings^[26]. Li, et al^[27] found that the thermal stress in segmented CoSb₃/Bi₂Te₃ device could be reduced by introducing a graded layer between CoSb₃ and copper electrode. However, all above researches treated multilayer systems as inert (without compositional or structural change) at high temperature. Actually, in TE joints, element diffusion and reaction are inevitable during whole service life. Therefore, the interfacial compositions and microstructures would continuously evolve, which influence the interfacial mechanical behavior unexpectedly.

This study reports a numerical analysis model based on finite element simulation method to investigate the dynamic interfacial stress at the SKD/DBL joint by taking the microstructure evolution into account. A single- layer model was established, and the interfacial stress in the TE joints with different DBL was calculated. Based on the experimental results on the interface observation and properties measuring of reaction layers, a multilayer model considering microstructure evolution was built to quantitively simulate the stress state of the aged TE joints. The tensile test results of SKD/Nb joints matched well with the simulation results.

A transient thermal-structural model is applied for stress in SKD/DBL joints to simulate the sintering- cooling process. Residual stress results from change of temperature and difference of material properties. For cooling period, the transient thermal conducting equation can be derived according to energy conservation law:

Where T is absolute temperature, ρ is the density, C_p is isobaric heat capacity, q is the heat flux vector, $\dot{q}$ is heat generation rate, [κ] stands for thermal conductivity matrix dependent with T. From Eq. (1-2), the temperature distribution in model can be worked out.

Once temperature decreases, the difference between CTE of SKD and Nb causes displacement, which can be converted to strain, and then results in stress according to generalized Hooke’s law. Thus the thermal-structural governing equation can be expressed as^[28]:

Where σ and ε are stress and strain vectors respectively, [D] is stiffness matrix, α is CTE vector, ΔT is the difference between present temperature and reference temperature (T₀). Finally, the stress distribution can be figured out.

For a certain site in materials, stress tensor can be divided to three principle stresses: σ₁, σ₂, σ₃, from large to small, respectively. According to the first strength theory, the maximum principle stress is the main reason for fracture, which coincides well with the fracture of brittle materials’ uniaxial tension. Therefore, the first principle stress σ₁ is chosen to evaluate stress intensity.

In initial state (t=0), adjacent materials are already bonded closely to form a zero-stress state. Thus, the initial temperature (960 K) is reference temperature (T₀) of thermal expansion. Then in cooling process, joint’s diameter is assumed to be unchanged, meanwhile upper and lower end could move vertically. Besides, model’s round side is considered as thermal insulated. On the end side surfaces, effect of cooling is also simplified as natural convection in air, in which the heat transfer coefficient was set to be 10 W·m²·K^-1.

SKD/Nb joints were fabricated and aged for tensile test. Yb_0.3Co₄Sb₁₂ (SKD) powders and Nb foil ((0.025± 0.015) mm) were loaded into a graphite die with a diameter of 50 mm, and then sintered by hot pressed for 90 min at 690 ℃ and 60 MPa under Ar atmosphere. Nb foil was placed between two SKD layers to form a sandwich structure. The obtained joint (ϕ50 mm×4 mm) was then cut into small cylinders (ϕ10 mm) and sealed in quartz ampules under vacuum. The sealed ampules were aged under 600 and 650 ℃ in furnace for various time, and the aged joints were denoted as “temperature-aging time”, such as 600-10 d. The tensile strength was measured by Instron-5566 universal testing system at room temperature. The microstructures of interface and fracture surface were observed by SEM (ZEISS Spura 55). The constituents of interface were measured by EDS (Oxford Instrument).

Nb, Mo, Zr, Ti were chosen as the DBL candidates for SKD^[22,29]. The effect of materials properties on interfacial mechanical stability was studied by using single- layer model. The CTE and Young’s modulus data obtained from COMSOL database are presented in Fig. 1(a-b). The calculated σ₁ are shown in Fig. 1(c). Compared with other DBL materials, the stress intensity of SKD/Mo is extremely high, which probably results from great difference of CTE between Mo and SKD (CTE=10.6× 10^-6 K^-1[20]). The stress of SKD/Zr interface is the lowest (1.1 GPa), even though Zr has the larger difference of CTE with SKD than that of Ti. Fig. 1(c) shows that the variation trend of interfacial stress is the same as that of Young’s modulus. It is observed that both CTE and Young’s modulus influence stress intensity, and large Young’s modulus induces high intensity stress. Moreover, in single-layer model, the thickness of DBL affects interfacial stress weakly. In consideration of reaction activity at high temperature^[22], Nb and Zr are employed for further research.

Actually, all the SKD joints undergo the elemental diffusion and chemical reaction during long-time service at high temperature^[8], which makes single-layer model inaccurate to describe the stress state of aged joints. SEM and EDS results of SKD/Nb joints after different aging time are listed in Fig. A2. No evident microcracks or micropores are observed in as-prepared joint. After aged at 600 ℃ for 5 d, NbSb₂ and CoSb₂ are detected, and micropores appear in the CoSb₂ layer. The thickness of NbSb₂ slightly increase with the increase of aging time (Fig. A2(b-d)), indicating that elevating temperature and prolonging aging time significantly aggravate diffusion and reaction. The reaction process can be described as following. At the initial stage of aging, SKD didn’t react with DBL, and there are only three layers, as the single-layer model is shown in Fig. 2(a). With aging accelerating, SKD reacted with Nb to form CoSb₂ and NbSb₂, and micropores appeared simultaneously because of the volume change. The reaction equation can be expressed in Eq. (4):

Thus, multilayer model was built to find out the influence of micropores and diffusion layers.

To simplify the modeling and calculation, there’re some requisite assumptions in the model: (1) CoSb₂ and NbSb₂ layers are assumed to be totally flat; (2) Micropores are of ellipsoid shape, in which x and y semi-major axis (a and b) are equal to a multiple of z semi-major axis (c); (3) Micropores locate on the interface between NbSb₂ and CoSb₂ periodically; (4) Diameter of cylindrical model is reduced; (5) Materials are seen as isotropic and completely linear-elastic; (6) Part of material properties (Table A1) is treated as constant value or simple function of temperature because of lack in experimental data. Combining materials’ molar mass and density listed in Table A1, relationships of thickness of each layer can be obtained, as shown in Table 1. Thus, total volume of micropores ought to equal the volume difference between models before and after aging. Therefore, the relationship (Eq. (5)) between average micropore size and thickness of NbSb₂ is listed:

Where r is radius of model, d_NbSb2 is thickness of NbSb₂, ΔV is total volume of micropores, m is the ratio of a (or b) to c, n is number of micropores at given total micropore volume. As long as micropores’ shape, numbers and positions are certain, d_NbSb2 is expected to decide the extent of aging. Micropores’ positions are considered to uniformly distribute at the CoSb₂/NbSb₂ interface. The number of micropores n was set as 1, 3, 7, and the distributions for each number are shown in Fig. 2(b-d). According to calculation results presented in Fig. A3, the ratio of a (or b) to c (m) made little difference to stress distribution, which won’t be discussed in the next part.

The diffusion and reaction process of SKD/Zr joint is similar with SKD/Nb, and the reaction layer is ZrSb₂^[22]. It can also be analyzed by using multilayer model, and the results are shown in Fig. A4.

The evolution of the interface can be simplified descripted by 3 parameters of n, d_DBL and d_NbSb2. The positions of micropores need to be set manually for each different n. When n changes from 1 to 7, the micropore size decreases from 5.3 μm to 1.2 μm. Furthermore, to eliminate the influence of abnormal stress caused by low-quality elements mesh in much thinner barrier layer and diffusion layer, average first principle stress is applied to evaluate interfacial stress intensity. Fig. 2(a) gives the initial stress distribution of SKD/Nb interface which is the same with the single layer model. The stress concentrates on Nb layer, and the maximum principle stresses of Nb/SKD interface is 2.383 and 1.46 GPa, respectively. Fig. 2(b-d) shows the σ₁ distribution in multilayer model with different micropore numbers n at the same given ΔV. After aging, all stress values with increasing n are much higher than the initial ones. For all cases, maximum stresses are found at the CoSb₂ layers. When d_NbSb2=10 μm, the maximum σ₁ are 5.83 GPa for n=1 and 7.79 GPa for n=7, indicating that the diffusion reaction induce large internal stress. Fig. 3 shows the d_DBL and d_NbSb2 dependent σ₁ at SKD/CoSb₂ and CoSb₂/NbSb₂ interfaces. The average stress of CoSb₂/NbSb₂ interface is always the largest value for different n, indicating that CoSb₂/NbSb₂ is the most unstable interface. In Fig. 3(b, e, h), contour lines are approximately parallel to d_DBL axis, which means the thickness of reactive NbSb₂ plays a dominate role on stress of CoSb₂/NbSb₂ interface.

SKD/Nb joints are fabricated and aged for tensile test. After fracture, bonding strength are calculated by Eq. (6):

Where σ_t is tensile strength, F_max is the maximum load, A is base area of joints. The results are listed in Table 2. It is obvious that accelerated aging worsens interfacial bonding severely. With the d_NbSb2 grows from 0 to 12 μm, the tensile strength decreases from 9.68 MPa to 1.46 MPa. The decreasing trend of tensile strength is consistent with the increasing trend of calculated stress. Furthermore, after tension, by comparing structures and compositions of the fracture surfaces, location of the weakest interface can be found. All of the aged joints break at CoSb₂ layer (Fig. 4(d-e)), while the unaged joints break at SKD/Nb interface (Fig. 4(a)). With elevated aging temperature or prolonged aging time, the proportions of CoSb₂ on fracture surface (SKD side, Fig. 4(d)) increases obviously (Table 2), indicating that the fracture locations tend to be the CoSb₂/NbSb₂ interface. The fracture locations are completely consistent with calculation results. All of the experimental results verify the validity of above simulation model.

In this study, a single-layer model was established to calculate the interfacial stress in the TE joints with different DBL. It’s found that the DBL materials with small modulus and similar CTE with SKD can reduce the interfacial stress. Based on the experimental results, a multilayer model considering microstructure evolution is built to quantitively simulate the stress state of the aged TE joints. Large thickness of reaction layers and volume changes can intensify stress at interface remarkably. Both in SKD/Zr and SKD/Nb joints, the biggest first principle stress locates at CoSb₂ layer. Tensile test results of SKD/ Nb joints fit simulation results well, proving the feasibility of this model to simulate the stress state in multilayer system containing complex microstructures, which is helpful to design the high stability electrode interface structure for SKD TE devices.

Supporting materials related to this article can be found at https://doi.org/10.15541/jim20190112.

[1] FAILAMANI F, BROZ P, MACCIÒ D, et al. Constitution of the systems {V,Nb,Ta}-Sb and physical properties of di-antimonides {V,Nb,Ta}Sb₂. Intermetallics, 2015, 65: 94-110.

[2] TAVASSOLI A, GRYTSIV A, FAILAMANI F, et al. Constitution of the binary M-Sb systems (M=Ti, Zr, Hf) and physical properties of MSb₂. Intermetallics, 2018, 94: 119-132.

[3] SALVADOR J R, YANG J, SHI X, et al. Transport and mechanical properties of Yb-filled skutterudites. Philosophical Magazine, 2009, 89(19): 1517-1534.

[4] GOTO Y, MIYAO S, KAMIHARA Y, et al. Electrical/thermal transport and electronic structure of the binary cobalt pnictides CoPn₂ (Pn=As and Sb). AIP Advances, 2015, 5(6): 067147.

[5] ZHAO D, LI X, JIANG W, et al. Fabrication of CoSb₃/MoCu thermoelectric joint by one-step SPS and evaluation. Journal of Inorganic Materials, 2009, 24(3): 545-548.

[6] BÖRNSTEIN L. CoSb₂: Crystal Structure, Physical Properties, in: Madelung U R O, Schulz M (Ed.), Non-tetrahedrally Bonded Binary Compounds ii. Berlin: Springer-Verlag, 2000.