Calculation of Thermal Expansion Coefficient of Rare Earth Zirconate System at High Temperature by First Principles

Abstract

Compounds of rare earth zirconates with pyrochlore structure are candidates for the application of thermal barrier coatings of next generation. In order to modify the mechanic properties and maintain the low thermal conductivity, other trivalent rare-earth element substitution is commonly used. Presently, investigation on the evaluation of the property of thermal expansion is attracting more attention. In this paper, a feature parameter of thermal expansion coefficient at high temperature (α_∞) was proposed by combining Grüneisen’s equation and the Debye heat capacity model. Using α_∞ model, the thermal expansion property of different compounds can be easily figured out by first principles. Firstly, α_∞ of ZrO₂, HfO₂, were calculated, and results are in good agreement with the experimental data from the literature. Moreover, α_∞ of La₂Zr₂O₇, Pr₂Zr₂O₇, Gd₂Zr₂O₇, and Dy₂Zr₂O₇ were calculated, and results demonstrated that the model of α_∞ is a useful tool to predict the thermal expansion coefficient at high temperature. Finally, Gd₂Zr₂O₇ with 4 different Yb dopant concentrations (Gd_1-xYb_x)₂Zr₂O₇ (x = 0, 0.125, 0.3125, 0.5) were calculated. Comparing with the experimental data from the literature, the calculation results showed the same tendency with the increasing of Yb concentration.

1. Introduction

It is well known that yttria-stabilized zirconia with a mass fraction of 8%(8YSZ) is widely used as top coating of thermal barrier coatings (TBCs) for aero engines [1,2]. However, with the continuously increasing demand of the thrust-to-weight ratio, 8YSZ is not workable because of its phase change and sintering since the temperature of the turbine front inlet is much greater than 1200 °C [3,4]. Gd₂Zr₂O₇, as the representative of rare earth zirconates, is evidently one of the most promising candidates for the application of next generation TBCs due to its lower thermal conductivity and higher phase stability [5,6]. However, it still suffers from the problem that mechanical properties are not high enough and thermal cycling performance is poor [7]. The thermal expansion properties play the key role. Comparatively, the thermal expansion coefficient of rare earth zirconates is about 9–10 × 10⁻⁶ K⁻¹(1073 K) [8], which is much lower than that of 8YSZ, about 11 × 10⁻⁶ K⁻¹(1073 K) [3]. The thermal expansion coefficient of NiCoCrAlY bonding layer, is about 17.5 × 10⁻⁶ K⁻¹(1273 K) [5]. The mismatch of the thermal expansion between the ceramic top layer and the bonding layer causes thermal stresses during thermal cycling, which can lead to cracks and failure of the TBCs system [9,10].

In order to improve the thermal expansion property of rare earth zirconates, doping with other trivalent rare earth elements is typically used. Excellent thermophysical properties such as high thermal stability, lower thermal conductivity, and high thermal expansion have been demonstrated in doped rare earth zirconates, such as Sm₂(Zr_1-xCe_x)₂O₇, (Nd_1-xGd_x)₂Zr₂O₇, Gd₂(Zr_1-xTi_x)₂O₇ and (Sm_1-xGd_x)Zr₂O₇ [6,11,12,13], etc.

Nevertheless, the research described above was carried out mainly by experiments, which are time costly and not fully adequate. For instance, heat conduction and thermal expansion of materials are closely related to lattice vibrations. However, to investigate the lattice vibrations of materials experimentally, neutron scattering or Raman spectroscopy is required [14]. Correspondingly, computational material simulations are more efficiently implemented to accelerate the design of materials. As we know, the first principles exhibit powerful capabilities in material design. Zhao [15] investigated the structure, mechanical properties, minimum thermal conductivity, and electronic properties of a series of Gd-site and Zr-site substituted Gd₂Zr₂O₇ pyrochlores by first principles. Atsushi Togo [16] calculated the thermal expansion properties of Ti₃SiC₂, Ti₃AlC₂, and Ti₃GeC₂ by the first principles combing with quasi-harmonic approximation (QHA). Feng [17] investigated the thermal expansion properties of rare earth zirconates (Ln₂Zr₂O₇, Ln = La, Nd, Sm and Gd) of pyrochlore structures by using first principles. Meanwhile, the calculation of properties of doped compounds by first principles is very difficult especially for the calculation of the thermal expansion property. For instance, the conventional cell of Gd₂Zr₂O₇ contains 88 atoms, and cell expansion must be performed to obtain solid solution structures with different doping concentrations, which makes the calculation of phonons extremely difficult.

Classically, Grüneisen proposed a thermal expansion equation [18], in which the thermal expansion coefficient(α) is determined by elastic properties and heat capacity (C_V). Further, C_V is the function of Debye temperature and temperature. When the temperature is much greater than Debye temperature, C_V can be considered as a constant. In this case, the calculation of α is simplified to the calculation of elastic properties, which makes the calculation much easier and faster. Coincidentally, the TBCs for aero engines operate at a high temperature, which is much higher than Debye temperature. Moreover, investigations showed that the coefficient of thermal expansion gradually increases with temperature increasing at high temperature, but the increasing rate gradually decreases. Reasonably, the thermal expansion coefficient at super high temperature(α_∞) can be used as a comparison factor to characterize the thermal expansion property of different dopant/concentration compound materials.

Therefore, in this paper, based on the Grüneisen’s equation, a computational model of α_∞ was developed, by which the calculation of a series of rare earth zirconates were implemented by the first principles.

2. Methodology

The supercell containing 88 atoms of Gd₂Zr₂O₇ was used for calculation. The structures of Yb doped Gd₂Zr₂O₇ were formed by replacing Gd atoms with different amounts of Yb atoms. The structural models for (Gd_1-xYb_x)₂Zr₂O₇ were built using the cluster expansion approach by calculating the lowest forming energy [19,20,21]. Further, the structures were optimized by the Birch–Murnaghan equation of state [22]. The elastic constants of the material were calculated by the stress–strain method [23].

The first principles calculations were based on density functional theory using the Vienna Ab initio Simulation Package (VASP) [24] with the generalized gradient approximation (GGA) for exchange-correlation energy, in the form of Perdew–Burke–Ernzerhof (PBE) [25]. The kinetic cut-off energy for the plane wave expansion was taken to be 600 eV in the Brillouin zone integrations using 2 × 2 × 2 k points. The average force acting on ions was reduced to 0.05 eV/Ang. Valence electrons included for distinct atoms were O 2s²2p⁴, Zr 5s¹4d³, Hf 5d³6s¹, Gd 6s²5p⁶5d¹, Yb 6s²5p⁶.

3. Results

3.1. Yb Doped Gd₂Zr₂O₇ Structure

The supercell of Gd₂Zr₂O₇ contains 88 atoms, including 16 Gd atoms, 16 Zr atoms and 56 oxygen atoms. A total of 2, 5 and 8 Yb atoms were used to replace the Gd atoms to obtain three different concentrations: (Gd_0.875Yb_0.125)₂Zr₂O₇, (Gd_0.875Yb_0.3125)₂Zr₂O₇, and (Gd_0.5Yb_0.5)₂Zr₂O₇, respectively. Correspondingly, the possible numbers of (Gd_1-xYb_x)₂Zr₂O₇ doped structures were $C_{16}^2$, $C_{16}^5$, and $C_{16}^8$. Excluding the equivalent structures, the number of unequal possible doped structures are 3, 35, and 97, separately. The forming energy E of each structure was calculated using the cluster expansion approach, according to Equation (1) [21], and the final doped structure was identified by the structure with the lowest forming energy.

$$E=\frac{\left(E_0-n_1\times E_1-n_2\times E_2\right)}{\left(n_1+n_2\right)}$$

(1)

wherein E₀ is the energy of the doped structure; E₁ and E₂ are the energy of the single cell of Gd₂Zr₂O₇ and Yb₂Zr₂O₇, and n₁ and n₂ are the number of dopant atoms, respectively. The calculation results were shown in Figure 1. According to the lowest forming energy, the geometrical configurations of three different doping structures were elaborated in Figure 2.

3.2. Lattice Constant and Elastic Modulus

Based on the geometrical configurations shown in Figure 2, the lattice constants and elastic modulus were calculated and listed in

Table 1. Lattice constant, the elastic constants (C₁₁, C₁₂, and C₄₄), bulk modulus (B), shear modulus (G), and Poisson’s ratio (μ) of rare earth zirconates.

	a0/(nm)	C11/(GPa)	C12/(GPa)	C44/(GPa)	B/(GPa)	G/(GPa)	G/B	μ
Gd2Zr2O7, cal.	1.056	316.4	106.7	84.2	176.6	91.9	0.52	0.278
Gd2Zr2O7, exp. [26]	1.054				174	93
(Gd0.875Yb0.125)2Zr2O7	1.055	312.3	100.2	83.4	170.6	91.7	0.54	0.272
(Gd0.6875Yb0.3125)2Zr2O7	1.052	308.4	96.7	83	167.2	91.6	0.55	0.269
(Gd0.5Yb0.5)2Zr2O7	1.050	310.4	96	82.5	167.7	91.7	0.55	0.269

.

It can be seen that the lattice constants decrease with the increasing of Yb dopant. The calculated bulk and shear modulus of Gd₂Zr₂O₇ are 176.6 Gpa and 91.9 GPa, respectively, which meets agreement with the data measured by experiments [26]. Generally, the bulk modulus and shear modulus of (Gd_1-xYb_x)₂Zr₂O₇ decrease with the increasing of Yb content, which are possibly caused by structure change. Subramanian M A [27] pointed out that doping of Yb atoms reduce the average cation radius ratio r(A³⁺)/r(B⁴⁺) and change the crystal structure of Gd₂Zr₂O₇ from pyrochlore to disorders in the structures. For the cubic phase, there are three independent elastic constants, C₁₁, C₁₂, and C₄₄ [28]. The calculated elastic constants are all positive, satisfying the generalized elastic stability criterion, namely, C₁₁ + 2C₁₂ > 0; C₄₄ > 0; C₁₁ − C₁₂ > 0, indicating that all studied structures are mechanically stable [29]. According to Pugh’s theory [30], when G/B < 0.5, the material is ductile; otherwise, the material is brittle. The G/B value of all materials calculated is greater than 0.5, indicating that they are brittle materials.

3.3. Thermal Expansion of Rare Earth Zirconates System

According to Grüneisen’s equation [18], the volumetric coefficient of thermal expansion (β) can be expressed as Equation (2).

$$\beta =\frac{\gamma }{B}\frac{C_V}{V}$$

(2)

Here, γ is Grüneisen’s constant, C_V is the heat capacity, B is the bulk modulus, and V is the molar volume.

As we know, the rare earth zirconate is in cubic phase. For cubic phase, γ is defined as function of the Poisson’s ratio (μ) as the following Equation (3) [31].

$$\gamma =\frac{3}{2}\left(\frac{1+\mu }{2-3\mu }\right)$$

(3)

C_V can be calculated by the Debye heat capacity model as the following Equation (4) [32]:

$$C_V=9N{k}_B{\left(\frac{T}{T_D}\right)}^3{\displaystyle {\int }_0^{\frac{T_D}{T}}\frac{x^4{e}^x}{{\left(e^x-1\right)}^2}}dx$$

(4)

wherein N = nN_A, N_A is Avogadro’s constant, and n is the number of atoms in the molecular formula; k_B is the Boltzmann constant; T_D is the Debye temperature. Meanwhile, T_D can be calculated by the following Equation (5) [33].

$$T_D=\frac{h}{k_B}{\left[\frac{3n}{4\pi }\left(\frac{N_A\rho }{M}\right)\right]}^{1/3}{\upsilon }_m$$

(5)

wherein h, ρ, M, υ_m is Planck’s constant, theoretical density, relative molecular mass, and speed of sound, respectively. υ_m is defined as Equation (6).

$${\upsilon }_m=\sqrt[3]{3}{\left(\frac{1}{{\upsilon }_L^3}+\frac{2}{{\upsilon }_S^3}\right)}^{-1/3}$$

(6)

υ_L, υ_S is the longitudinal and transverse sound velocity and can be expressed as Equations (7) and (8), separately.

$${\upsilon }_L=\sqrt{\frac{3B+4G}{3\rho }}$$

(7)

$${\upsilon }_S=\sqrt{\frac{G}{\rho }}$$

(8)

For cubic crystals, β = 3α, α is the linear expansion coefficient [17]. From the above Equations (4)–(8), it can be deduced that α can be expressed as Equation (9):

$$\alpha =\frac{9N{k}_B\frac{1+\mu }{2-3\mu }{\left(\frac{T}{T_D}\right)}^3{\displaystyle {\int }_0^{\frac{T_D}{T}}\frac{x^4{e}^x}{{\left(e^x-1\right)}^2}}dx}{2BV}$$

(9)

Moreover, Equation (9) can be rewritten as Equation (10).

$$\alpha =\left(\frac{9N{k}_B\frac{1+\mu }{2-3\mu }{T}^3{\displaystyle {\int }_0^{\frac{T_D}{T}}\frac{x^4{e}^x}{{\left(e^x-1\right)}^2}}dx}{2BV}\right)T_D^{-3}$$

(10)

According to the calculation result shown in Table 1, B, μ, and V of different (Gd_1-xYb_x)₂Zr₂O₇ vary within 4% difference, which can be considered as constant. Thus, it is indicated from Equation (10) that α is proportional to T_D⁻³ at the same temperature, which is compliant with Ruffa’s equation [34].

Furthermore, supposing that the temperature is much greater than the Debye temperature, the integral term in Equation (9) can be mathematically simplified and finally Equation (11) can be obtained.

$${\alpha }_{\infty }=\frac{3N{k}_B\frac{1+\mu }{2-3\mu }}{2BV}$$

(11)

Here, α_∞ representatives the linear thermal expansion coefficient at super high temperature. Actually, the TBCs are working under the temperature (e.g., the temperature of combustion chamber in F135 turbine engine can be up to 2253 K [35]) much higher than T_D of rare earth zirconate (about 500 K) [36]. Another one, α of rare earth zirconate, increases with the increase in temperature, meanwhile the increasing rate slows down more and more [37,38,39]. Therefore, α_∞ can be likely used to compare the difference of thermal expansion property of different dopant/concentration for the same compound.

3.4. The Validity of α_∞ Model

Cubic ZrO₂ and HfO₂ are of typical fluorite structure which is the same as that of rare earth zirconates [40]. Firstly, the lattice parameter and elastic properties were calculated by first principles. Secondly utilizing the α_∞ model, the thermal expansion properties of cubic ZrO₂ and HfO₂ were calculated, and both results were listed in

Table 2. Lattice constant, bulk modulus (B), shear modulus (G), Poisson’s ratio (μ), and α_∞ of ZrO₂ and HfO₂.

	a0/(nm)	B/(GPa)	G/(GPa)	μ	α∞ (K−1)
ZrO2, cal.	0.512	238.5	100.6	0.316	9.72 × 10−6
ZrO2 [a]	0.515	235	103	0.31
HfO2, cal.	0.508	253.8	112.6	0.307	9.05 × 10−6
HfO2 [b]	0.508	248	115	0.3

[a] Materials data on ZrO₂ (SG:225) by Materials Project. ID:mp-1565. [b] Materials data on HfO₂ (SG:225) by Materials Project. ID:mp-550893.

. Comparing to the data from the material project database, calculation results of lattice parameter and elastic properties are very close to the same level. Hong [41] and Irshad, K.A. [42] measured the linear thermal expansion coefficient (α) of cubic ZrO₂ and HfO₂ by in situ high-temperature X-ray diffraction, and it is (12 ± 3) × 10⁻⁶ K⁻¹ and 8.80 × 10⁻⁶ K⁻¹, respectively. Comparably, the calculated results of α_∞ are 9.72 × 10⁻⁶ K⁻¹ and 9.05 × 10⁻⁶ K⁻¹. It is revealed that both are a good match.

In order to further verify the validity of the α_∞ model, α_∞ of serial rare earth zirconates, including La₂Zr₂O₇, Pr₂Zr₂O₇, Gd₂Zr₂O₇, and Dy₂Zr₂O₇ were calculated. The data of elastic property was cited from the literature [36], and both were shown in

Table 3. Lattice constant a₀, bulk modulus (B), shear modulus (G), Poisson’s ratio (μ) of rare earth zirconates system [36].

	a0/(nm)	B/(GPa)	G/(GPa)	μ	Thermal Expansion Coefficient/(10−6K−1)
					α∞	α/(1073 K) [8]
La2Zr2O7	1.081	176	87	0.302	9.755	8.883
Pr2Zr2O7	1.072	155	103	0.26	9.857	9.415
Gd2Zr2O7	1.052	165	63	0.284	10.61	10.094
Dy2Zr2O7	1.054	164	90	0.268	10.057	9.166

. α_∞ and α were plotted in Figure 3. Compared to α measured by experiment [8], α_∞ are clearly higher because α were measured at 800 °C which is not too much higher than Debye temperature. Henry Lehmann [37] measured the thermal expansion coefficients of Gd₂Zr₂O₇ and La₂Zr₂O₇, which were 10.652 × 10⁻⁶ K⁻¹ (1473 K) and 9.09 × 10⁻⁶ K⁻¹ (1373 K), respectively. The results are very close to α_∞ calculated. It is further demonstrated that α_∞ can be a useful tool to predict the thermal expansion coefficient at high temperature.

3.5. The Effect of Yb Doping of Gd₂Zr₂O₇ on α_∞

α_∞ of Gd₂Zr₂O₇ with 4 different Yb doping contents Gd₂Zr₂O₇, (Gd_0.875Yb_0.125)₂Zr₂O₇, (Gd_0.875Yb_0.3125)₂Zr₂O₇, and (Gd_0.5Yb_0.5)₂Zr₂O₇ were calculated.

Table 4. Density (ρ), longitudinal wave velocity (υ_L), shear wave velocity of sound (υ_S), average velocity of sound (υ_m) and Debye temperature (T_D) of rare earth zirconates.

	ρ/(kg·m−3)	υL/(m·s−1)	υS/(m·s−1)	υm/(m·s−1)	TD/(K)
Gd2Zr2O7	6868	6600	3659	4075	511
(Gd0.875Yb0.125)2Zr2O7	6944	6496	3635	4046	508
(Gd0.6875Yb0.3125)2Zr2O7	7059	6402	3602	4007	504
(Gd0.5Yb0.5)2Zr2O7	7176	6357	3574	3977	502

shows the calculation result of theoretical density, sound velocity, and Debye temperature.

It can be seen that T_D of Gd₂Zr₂O₇ was calculated to be 511 K, which is in good agreement with that measured by Toshiaki Kawano [8]. T_D is decreases with the increase of Yb dopant, which is dominated by the decrease of average velocity of sound (υ_m).

Figure 4 shows the difference between calculated α_∞ and α measured by experiment at 1073 K [43].

It is clear that Yb doping can increase the thermal coefficient greatly because theoretical density increases meanwhile T_D decreases with Yb doping, as shown in Table 4. However, α_∞ is a bit lower than α. Feng [17] also discussed and explained the problem. Actually, the first principles were developed based on the material being from an ideally perfect crystal. However, for real bulk materials, various defects (e.g., vacancies and dislocations) and pores ineluctably existed. In general, the total energy of a crystal with defects is higher than that of an ideally perfect crystal, and the anharmonic effect may be affected by various defects and pores in the structure. In addition, the density of the tested ceramic coupon is certainly lower than that of theoretical density. That is the reason why the thermal expansion coefficient measured by experiment is higher than that calculated.

Figure 4 also reveals that both α_∞ and α increase with an increase in Yb content. The increase of thermal coefficient is more remarkable with lower Yb doping concentration. With higher Yb doping concentration, the growth rate slows down. Both the measured and calculated results have the same tendency.

Thermal expansion is related to the crystal structure and the electronic structure [44]. The PDOS of (Gd_1-xYb_x)₂Zr₂O₇ crystals are shown in Figure 5. In terms of bonding, Gd/Yb-5d, Zr-4d, and O-2p overlap, which means that the electronic state around the Fermi level is primarily determined by the relatively weak p-d bond between O-2p and Zr-4d (or Gd/Yb-5d). The main change from doping is the relative position changes in Gd/Yb- 5d, which may be caused by the difference in the valence electron layers between Gd and Yb. With the increase of Yb content, the total state density curve moves slightly to the lower energy level. Low crystal energy means high coefficient of thermal expansion [45], and the rare earth element Yb influences the O-2p states due to hybridization [44]. The p-d bond strength decreased with the increase in Yb content. According to the value of the PDOS ordinate, the PDOS of Zr-4d is the largest, so the p-d bond strength of Zr-O is higher than Yb-O and Gd-O.

4. Conclusions

Rare earth zirconates are candidates for next generation TBCs, and it is important to develop a method to characterize the thermal expansion property. Combining the Grüneisen’s equation and the Debye heat capacity model, an efficient model of α_∞ to characterize the thermal expansion coefficient at super high temperature was established. Firstly, using the α_∞ model, the high temperature thermal expansion coefficients of cubic ZrO₂ and cubic HfO₂ were calculated to be 9.72 × 10⁻⁶ K⁻¹and 9.05 × 10⁻⁶ K⁻¹, respectively, which are in agreement with those shown in the literature. Secondly, α_∞ of serial rare earth zirconates, including La₂Zr₂O₇, Pr₂Zr₂O₇, Gd₂Zr₂O₇, and Dy₂Zr₂O₇ were calculated, and results demonstrated that α_∞ can be a useful tool to predict the thermal expansion coefficient at high temperature. Lastly, α_∞ of (Gd_1-xYb_x)₂Zr₂O₇ with four different doping contents were calculated, and results showed the same tendency as that measured by experiments. Generally, by characterizing the thermal expansion coefficient at high temperature through the elastic properties and Debye temperature of the material, the complicated calculation of phonon spectrum can be avoided. Thus, the model of α_∞ has the broad application prospect to predict the thermal expansion property at high temperature for other rare earth zirconates.