1 Introduction

Due to the special ability in converting energy between heat and electricity, thermoelectric materials (TEMs) have many engineering applications, such as waste heat recovery^[1], thermoelectric (TE) cooling^[2] and thermal protection system in the aerospace industry^[3]. However, the low conversion efficiency limits the development of TEMs. One available method to improve the conversion efficiency is to use multilayered TEMs^[4]. Therefore, the study of multilayered TEMs becomes a research interest.

Extensive researches about multilayered TEMs have been carried out. For example, buckling of TEMs with a large length-to-thickness ratio was analyzed^[5]. The buckling axial forces of thin-film TE generator made up of pn-junction were given^[6]. The problem of axial thermal stress^[7] and interlaminar thermal stress in layered TEMs^[8] was studied. The conversion efficiency is a key parameter for TEMs. Energy conversion efficiency for TE devices made up of pn-junction with temperature dependent and independent electrical conductivity was analyzed^[9]. Similarly, the efficiency of sandwiched TEMs with a graded interlayer whose properties vary with temperature was studied^[10]. By taking the Peltier effect into account, the overall macro-performance of multilayered TEMs was discussed^[11]. The effect of interface layers on the performance of annular TE generators was carried out^[12].

Most TEMs are brittle in nature. Therefore, understanding their fracture properties is of great importance and will contribute to the development of multilayered TEMs^[13-14]. However, because of the natures of non-homogeneity and poor interfacial shear strength, multilayered material systems are prone to damage by delamination^[15]. The delamination may result in the loss of functionality of TE device. Therefore, mechanics (especially the fracture mechanics) research of multilayered TEMs is very important. However, fracture mechanics research of thermally induced delamination of multilayered TEMs has not yet been conducted.

The fracture mechanics analysis of delamination in a TE pn-junction sandwiched by an insulating layer was conducted. Time-varying energy release rates at the left and right tips of the delamination crack are given. The critical temperature difference that causes the delamination propagation is identified. Effects of temperature difference, electric current, and thickness of the insulating layer on the energy release rate are discussed.

2 Statement of the problem

Shown in Fig. 1 is a three-layer TEM which consists of an n-type TE layer and a p-type TE layer separated by an insulating layer. In the n-type TEM, the electric current is propagated by electrons. The charge carriers in the p-type TEM are holes which travel in the opposite direction with electrons^[16]. The two TE layers have the same material properties and thickness h^[5]. The thickness of the insulating layer is denoted by H. A through delamination crack occurs between the n-type layer and the insulating layer at the region l < x < l+a, where a is the length of the delamination crack. Expect that, for the region l < x < l+a, all layers are perfectly bonded together. The delamination part is replaced by equivalent axial forces F_L and F_R. The subscripts L and R denote the regions -C < x < l and l+a < x < C, respectively. The axial force of the delamination part is assumed to be uniform, which is presented by the distributed load χ. The electric potentials of the TEM of Fig. 1 are V=V₀ at the left end x=-C and V=V_h at the right end x=C. The initial temperature of the entire system is T₀, and the temperature at the right end is suddenly changed to T_h (assume T_h > T₀ in this paper), where the subscripts h and 0 denote the hot end and the cold end of the TEM, respectively. This means that the left end and the right end of the TEM of Fig. 1 are subject to a combination of electric potential and temperature, respectively. Furthermore, the boundary conditions of the surfaces of delamination crack, the upper border , and the lower border of the TEM are regarded as electrical and thermal insulations.

Since the thickness of each layer is very small compared with its length, the heat flux and stress fields only vary along the length direction but are constant along the thickness and width directions of the system. Therefore, the problem is one-dimensional. The constitutive equations of the TEMs when there is no free electric charge and heat source are^[8]

(1)

(2)

where J and q are the electric current density and heat flux, respectively, σ, s, and k are, respectively, the electric conductivity, the Seebeck coefficient, and the heat conductivity of the TE layers in which i (i=N, P) indicate the n-type TE layer and the p-type TE layer, respectively. For the transient-state problem, the equilibrium equations are and , where ρ_i and t are the thermal diffusion coefficient of the TE layer and the time.

3 Axial force

The transient-state temperature field and axial force of multilayered TEMs without delamination crack were presented by Wang and Cui^[8]. A brief review of the temperature field and axial force is given as follows. By taking Fourier's heat conduction, Joule's heat, and the energy accumulation into account, the following temperature governing equation of the multilayered TEMs is given:

(3)

By considering the energy balance and thermal boundary conditions (T=T₀ at x=-C and T=T_h at x=C), the effective transient temperature field of the TEMs is given as

(4)

where , β_n=(2n-1)π/(2L), n=1, 2, ..., and ρ_eq and k_eq are the effective thermal diffusion coefficient and the effective thermal conductivity, respectively, defined as ρ_eq=ρ_N+ρ_IH/(2h) and k_eq=k_N+k_IH/(2h)^[10]. The subscript I denotes the insulating layer. ΔT is the temperature difference between the hot end and the cold end of the system, i.e., ΔT=T_h-T₀.

The transient axial force F_N(x, t) in the n-type layer is

(5)

where A=-a₂C/sinh(ηC), B₁=-(a₁C²+a₃)/cosh(ηC), B₂=(-1)ⁿCd₁/cosh(ηC), a₁=3γJ²/(2η²k_eqσ), a₂=-γΔT/(η²C), a₃=-γ(J²C²/(k_eqσ) +ΔT)/(2η²)+2a₁/η², d₄=(2β_nd₁-γD₁)/(η²+β_n²), d₁=γD₁β_n/(η²+β_n²), η²=(λ_N+2λ_I)/(κ_I+κ_N), and γ =(α_N-α_I)/(κ_I+κ_N). α_i (i=N, I) are the thermal expansion coefficients of the n-type layer and the insulating layer, respectively. λ_i and κ_i are the coefficients of axial and interfacial compliances, respectively, defined as λ_i=(1-υ_i)/(E_ih_i) and κ_i=2h_i(1-υ_i)/(3E_i). υ_i, E_i, and h_i are Poisson's ratio, Young's modulus, and the thickness, respectively.

Since the delamination crack is parallel to the direction of heat flux, the temperature field given by Eq. (4) will not be perturbed by the delamination crack. Therefore, we will use the temperature field given by Eq. (4) to study the transient axial forces F_N^*(x, t) and F_I^*(x, t). The superscript * denotes the state of the n-type layer delaminating from the insulating layer. Since the n-type layer and the p-type layer are symmetric with respect to x, the displacement and axial force in the two TE layers are the same. To obtain the solutions of transient axial force, we first present the axial displacements u_N^- (x) and u_I⁺ (x). The superscripts -- and + denote, respectively, the displacement at the lower extreme fiber of the n-type layer and that at the upper extreme fiber of the insulating layer. The expressions of u_N^- (x) and u_I⁺ (x) are given as follows^[8]:

(6a)

(6b)

where τ_NI^*(x, t) is the shear stress which acts at the interface between the n-type layer and the insulating layer. By making use of the force balance condition along the x-direction, one obtains 2F_N^*(x, t)+F_I^*(x, t)=0. Setting the displacement u_N^- (x) equal to u_I (x) yields

(7)

Using the equation and differentiating Eq. (7) with respect to x, we can get . Under the mechanical boundary conditions (x=± C, F_N^* (x, t)=0), the solution of F_N^* (x, t) can be presented as

(8)

where the subscript m=L, R, and and are the unknown constants that can be solved from the mechanical boundary conditions.

The equivalent axial forces F_L and F_R can be directly obtained by applying F_L =F_N (l, t) at x=l and F_R =F_N (l+a, t) at x=l+a in Eq. (5). Substitution of the mechanical boundary conditions x=-C, F_N^*(-C, t)=0 and x=l, F_N^*(l, t)=F_L into Eq. (8), respectively, the unknown coefficients and can be obtained. Similarly, under the mechanical boundary conditions x=l+a, F_N^*(l+a, t)=F_R and x=C, F_N^*(C, t)=0, the unknown constants and can also be obtained. To obtain the axial force of the delamination part, we need to calculate the distributed load q_u firstly. By considering the force balance condition of the delamination part (l < x < l+a), χ =(F_R -F_L)/a can be given. Therefore, the axial force of the delamination part is F_N^*(x, t)=F_L+χ(x-l). The axial force of the insulating layer can be directly given by F_I^*(x, t)=-2F_N^*(x, t).

4 Delamination energy release rate

The energy release rate can be used as a criterion to judge whether the delamination propagates or not. Since the applied load is zero, the total potential energy equals the elastic strain energy U. The energy release rate , where dS is the incremental fracture area, dS=2da. The elastic strain energy per unit length of Fig. 1 is

(9)

where A_N and A_I represent the cross areas of the n-type layer and the insulating layer, respectively. By using 2F_N^*(x)+F_I^*(x)=0, Eq. (9) can be rewritten as

(10)

Integrating Eq. (10) with respect to x at the interval [-C, l+a] and [l+a, C], then substituting it into the equation , respectively, we can get the following energy release rates at the left tip and the right tip of the delamination crack:

(11a)

(11b)

Equations (11a) and (11b) indicate that the energy release rates G_L and G_R are determined by the equivalent axial forces F_L and F_R. The equivalent axial forces can be evaluated by Eq. (5), in which the applied electric current density J is not coupled with the temperature difference ΔT across the system. Therefore, we can investigate the effect of J and ΔT on delamination energy release rates G_L and G_R separately.

5 Numerical results and discussion

For the purpose of numerical demonstration, the TE layer and the insulating layer are selected as the Bi₂Te₃ and the ceramic, respectively. The material properties used in this section are α_N =1.68× 10^-5 K, α_I =0.65× 10^-5 K, E_N =47 GPa, E_I =380 GPa, υ_N =0.4, υ_I=0.26, k_N =2.2 W/(m·K), k_I=17 W/(m·K), ρ_N =1.312× 10⁶, ρ_I =3.176× 10⁶, σ_N =1×10⁵ Ω^-1·m^-1, and s_N=200 μV/K^[8]. The semi-length of the structure is C=150 μm. Unless otherwise specified, the TE layer and the insulating layer are assumed to have the same thickness h=H=3 μm. As presented in Eq. (8), the solution of axial force is given by the sum of an infinite series. All results in this section are for the truncation number 50, which has resulted in the convergent results. The delamination crack is assumed to be symmetric .

5.1 Model verification

In this section, the numerical model is verified by the finite element software ANSYS. The temperatures at the cold end and the hot end are, respectively, T₀=290 K and T_h=330 K. The electric potentials are prescribed as V₀=0 at x=-C and V_h=0.01 V at x=C. For the one-dimensional problem, the electric current density J is a constant, which can be expressed as J=-σ (V_h-V₀+sΔT)/(2C)=6× 10⁶ A/m^2[8]. For the typical delamination length a=10 μm, Fig. 2 illustrates the axial stress in the delamination region of the multilayered TEMs. Figure 2 indicates that the axial force in the delamination part is uniform.

The crack tip intensity factor K at the left tip and the right tip of delamination crack is also analyzed in ANSYS. The relation between the crack intensity factor and the energy release rate for the case of plane stress is G=K²/E. For the typical length of delamination crack a=5 μm, 10 μm, and 30 μm, Table 1 lists the energy release rates from the numerical model and ANSYS. It can be seen that the numerical results agree well with the finite element results.

5.2 Effect of the temperature difference on the delamination energy release rate

In this section, the applied electric current density is prescribed as J=0. Substituting J=0 into Eqs. (11a) and (11b) gives the delamination energy release rates G_LT and G_RT. The subscript T denotes the effect of temperature difference. To normalize delamination energy release rates G_LT and G_RT, we introduce a reference energy release rate G_0T=(ΔTγ/η²)²(1/(E_Nh) +2/(E_IH)). Figure 3 presents the effect of temperature difference ΔT on the normalized G_LT and G_RT. The solid curves denote the distributions of G_LT. The dashed curves denote the distributions of G_RT. It is found that the value of the energy release rate at the right tip of the delamination is larger than that at the left tip. Overall, the values of maximum G_LT and G_RT at the transient state are bigger than those of the steady state (e.g., t≥0.1 s). This identifies the significance of considering the transient state for the analysis of delamination of multilayered TE structures.

5.3 Effect of the electric current on the delamination energy release rate

In this section, the temperature difference is prescribed as ΔT=0. With substitution of ΔT=0 into Eqs. (11a) and (11b), one can obtain the delamination energy release rates G_LJ and G_RJ. The subscript J denotes the effect of electric current. The reference energy release rate G_0J=(J²C²γ/(η²k_eqσ))²(1/(hE_N) +2/(HE_I)) is used to normalize G_LJ and G_RJ.

Figure 4 presents the effect of electric current on the delamination energy release rate. Again, the solid and dashed curves denote the energy release rates at the left and right tips of delamination crack, respectively. It is found that the value of maximum energy release rate decreases with time. The energy release rate at the right tip of delamination crack is larger than that at the left tip.

5.4 Effect of the insulating layer thickness on the delamination energy release rate

Equation (3) clearly indicates that, though the constitutive equation of Eq. (1) is one-dimensional, the ratio of the insulating layer thickness to the TE layer thickness has a great effect on the distribution of temperature. As described in Section 4, the temperature greatly contributes to the distribution of G. This means that the thickness of insulating layer can affect the distribution of G. The thermal and electric boundary conditions applied at the TE system are the same as those of Subsection 5.1. For the typical time t=0.1 s, Fig. 5 presents the normalized G_R for different values of insulating layer thickness H=30 μm, 3 μm, 0.6 μm, and 0.1 μm. The reference energy release rate is G₀=h/(1/E_N+2/E_I). It is found that a thinner insulating layer results in a smaller G_R. This means that the thickness of the insulating layer should be as thin as possible in engineering application.

5.5 Critical temperature difference for the delamination propagation

Since the delamination energy release rate at the right tip is bigger than that at the left tip, G_R will be used to investigate the critical delamination temperature difference ΔT_c. Based on the energy release rate criterion, the delamination crack will propagate when the expression G_R ≥ K_c²/E_I is satisfied, in which K_c is the interface fracture toughness between the p-type layer and the insulating layer. With the substitution of energy release rate given by Eq. (11b), one can obtain ΔT_c. In this section, the value of K_c is chosen as 0.06 MPa· m^1/2[17]. For the typical time t=0.1 s and the given fracture toughness K_c, Fig. 6 shows the critical delamination temperature difference for different thicknesses of the insulating layer. It can be seen that ΔT_c decreases with the length of the delamination crack. For the same length of delamination crack, a thinner insulating layer results in a higher ΔT_c.

6 Conclusions

The fracture mechanics analysis of delamination in TE systems of the n-type layer and p-type layer separated by an insulating layer is carried out. The time-varying delamination energy release rate is given. Effects of the temperature difference, the applied electric current, and the thickness of insulating layer on the delamination energy release rate are discussed. The critical temperature difference of delamination propagation is explored by using an energy release rate criterion. It is found that the maximum delamination energy release rate of the transient state is bigger than that of the steady state. This suggests that considering the transient state when constructing the stability design of TE systems is essential. The thinner insulating layer results in the smaller delamination energy release rate and the higher critical delamination temperature difference. Therefore, the thickness of insulating layer should be as thin as possible in engineering applications.