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

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 ZrO2, HfO2, were calculated, and results are in good agreement with the experimental data from the literature. Moreover, α∞ of La2Zr2O7, Pr2Zr2O7, Gd2Zr2O7, and Dy2Zr2O7 were calculated, and results demonstrated that the model of α∞ is a useful tool to predict the thermal expansion coefficient at high temperature. Finally, Gd2Zr2O7 with 4 different Yb dopant concentrations (Gd1-xYbx)2Zr2O7 (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]. Gd2Zr2O7, 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−6 K−1(1073 K) [8], which is much lower than that of 8YSZ, about 11 × 10−6 K−1(1073 K) [3]. The thermal expansion coefficient of NiCoCrAlY bonding layer, is about 17.5 × 10−6 K−1(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 Sm2(Zr1-xCex)2O7, (Nd1-xGdx)2Zr2O7, Gd2(Zr1-xTix)2O7 and (Sm1-xGdx)Zr2O7 [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 Gd2Zr2O7 pyrochlores by first principles. Atsushi Togo [16] calculated the thermal expansion properties of Ti3SiC2, Ti3AlC2, and Ti3GeC2 by the first principles combing with quasi-harmonic approximation (QHA). Feng [17] investigated the thermal expansion properties of rare earth zirconates (Ln2Zr2O7, 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 Gd2Zr2O7 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). Further, C is the function of Debye temperature and temperature. When the temperature is much greater than Debye temperature, C 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 Gd2Zr2O7 was used for calculation. The structures of Yb doped Gd2Zr2O7 were formed by replacing Gd atoms with different amounts of Yb atoms. The structural models for (Gd1-xYbx)2Zr2O7 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 2s22p4, Zr 5s14d3, Hf 5d36s1, Gd 6s25p65d1, Yb 6s25p6.

3. Results

3.1. Yb Doped Gd2Zr2O7 Structure

The supercell of Gd2Zr2O7 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: (Gd0.875Yb0.125)2Zr2O7, (Gd0.875Yb0.3125)2Zr2O7, and (Gd0.5Yb0.5)2Zr2O7, respectively. Correspondingly, the possible numbers of (Gd1-xYbx)2Zr2O7 doped structures were , , and . 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. wherein E0 is the energy of the doped structure; E1 and E2 are the energy of the single cell of Gd2Zr2O7 and Yb2Zr2O7, and n1 and n2 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.

Figure 1

Calculation of E of (Gd1-xYbx)2Zr2O7 (x = 0.125, 0.3125, 0.5).

Figure 2

Geometrical configurations of (Gd1-xYbx)2Zr2O7, (a) Gd2Zr2O7; (b) (Gd0.875Yb0.125)2Zr2O7; (c) (Gd0.6875Yb0.3125)2Zr2O7; (d) (Gd0.5Yb0.5)2Zr2O7.

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.

Table 1

Lattice constant, the elastic constants (C11, C12, and C44), 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 Gd2Zr2O7 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 (Gd1-xYbx)2Zr2O7 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(A3+)/r(B4+) and change the crystal structure of Gd2Zr2O7 from pyrochlore to disorders in the structures. For the cubic phase, there are three independent elastic constants, C11, C12, and C44 [28]. The calculated elastic constants are all positive, satisfying the generalized elastic stability criterion, namely, C11 + 2C12 > 0; C44 > 0; C11 − C12 > 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).

Here, γ is Grüneisen’s constant, C 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].

C can be calculated by the Debye heat capacity model as the following Equation (4) [32]: wherein N = nN, N is Avogadro’s constant, and n is the number of atoms in the molecular formula; k is the Boltzmann constant; T is the Debye temperature. Meanwhile, T can be calculated by the following Equation (5) [33]. wherein h, ρ, M, υ is Planck’s constant, theoretical density, relative molecular mass, and speed of sound, respectively. υ is defined as Equation (6).

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

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):

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

According to the calculation result shown in Table 1, B, μ, and V of different (Gd1-xYbx)2Zr2O7 vary within 4% difference, which can be considered as constant. Thus, it is indicated from Equation (10) that α is proportional to T 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.

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 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 ZrO2 and HfO2 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 ZrO2 and HfO2 were calculated, and both results were listed in Table 2. 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 ZrO2 and HfO2 by in situ high-temperature X-ray diffraction, and it is (12 ± 3) × 10−6 K−1 and 8.80 × 10−6 K−1, respectively. Comparably, the calculated results of α∞ are 9.72 × 10−6 K−1 and 9.05 × 10−6 K−1. It is revealed that both are a good match.

Table 2

Lattice constant, bulk modulus (B), shear modulus (G), Poisson’s ratio (μ), and α∞ of ZrO2 and HfO2.

	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 ZrO2 (SG:225) by Materials Project. ID:mp-1565. [b] Materials data on HfO2 (SG:225) by Materials Project. ID:mp-550893.

In order to further verify the validity of the α∞ model, α∞ of serial rare earth zirconates, including La2Zr2O7, Pr2Zr2O7, Gd2Zr2O7, and Dy2Zr2O7 were calculated. The data of elastic property was cited from the literature [36], and both were shown in Table 3. α∞ 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 Gd2Zr2O7 and La2Zr2O7, which were 10.652 × 10−6 K−1 (1473 K) and 9.09 × 10−6 K−1 (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.

Table 3

Lattice constant a0, 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

Figure 3

α∞ and α of La2Zr2O7, Pr2Zr2O7, Gd2Zr2O7, and Dy2Zr2O7.

3.5. The Effect of Yb Doping of Gd2Zr2O7 on α∞

α∞ of Gd2Zr2O7 with 4 different Yb doping contents Gd2Zr2O7, (Gd0.875Yb0.125)2Zr2O7, (Gd0.875Yb0.3125)2Zr2O7, and (Gd0.5Yb0.5)2Zr2O7 were calculated. Table 4 shows the calculation result of theoretical density, sound velocity, and Debye temperature.

Table 4

Density (ρ), longitudinal wave velocity (υ), shear wave velocity of sound (υ), average velocity of sound (υ) and Debye temperature (T) 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

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

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

Figure 4

α∞ and α of (Gd1-xYbx)2Zr2O7.

It is clear that Yb doping can increase the thermal coefficient greatly because theoretical density increases meanwhile T 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 (Gd1-xYbx)2Zr2O7 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.

Figure 5

Partial density of states of (Gd1-xYbx)2Zr2O7: (a) x = 0; (b) x = 0.125; (c) x = 0.3125; (d) x = 0.5.

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 ZrO2 and cubic HfO2 were calculated to be 9.72 × 10−6 K−1and 9.05 × 10−6 K−1, respectively, which are in agreement with those shown in the literature. Secondly, α∞ of serial rare earth zirconates, including La2Zr2O7, Pr2Zr2O7, Gd2Zr2O7, and Dy2Zr2O7 were calculated, and results demonstrated that α∞ can be a useful tool to predict the thermal expansion coefficient at high temperature. Lastly, α∞ of (Gd1-xYbx)2Zr2O7 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.