Intrinsic Thermal Shock Behavior of Common Rutile Oxides

Abstract

Rutile TiO₂, VO₂, CrO₂, MnO₂, NbO₂, RuO₂, RhO₂, TaO₂, OsO₂, IrO₂, SnO₂, PbO₂, SiO₂, and GeO₂ (space group P4₂/mnm) were explored for thermal shock resistance applications using density functional theory in conjunction with acoustic phonon models. Four relevant thermomechanical properties were calculated, namely thermal conductivity, Poisson’s ratio, the linear coefficient of thermal expansion, and elastic modulus. The thermal conductivity exhibited a parabolic relationship with the linear coefficient of thermal expansion and the extremes were delineated by SiO₂ (the smallest linear coefficient of thermal expansion and the largest thermal conductivity) and PbO₂ (vice versa). It is suggested that stronger bonding in SiO₂ than PbO₂ is responsible for such behavior. This also gave rise to the largest elastic modulus of SiO₂ in this group of rutile oxides. Finally, the intrinsic thermal shock resistance was the largest for SiO₂, exceeding some of the competitive phases such as Al₂O₃ and nanolaminated Ti₃SiC₂.

1. Introduction

Rutile oxides (space group P4₂/mnm, prototype TiO₂), including TiO₂, VO₂, CrO₂, MnO₂, NbO₂, RuO₂, RhO₂, TaO₂, OsO₂, IrO₂, SnO₂, PbO₂, SiO₂, and GeO₂, are very common oxides and broadly explored due to their interesting properties [1,2]. For example, TiO₂ possesses a large band gap of about 3 eV [3]. RuO₂ is electrically conductive [4], which is highly unusual for oxides. On the one hand, NbO₂ exhibits the highest known Mott transition temperature of approximately 800 °C [5,6,7]. On the other hand, VO₂ undergoes the Mott transition at a low temperature of 68 °C [8], which is relevant for some applications such as smart windows. SiO₂, in its various forms, is known for high thermal shock resistance [9]. In general, oxides are refractory solids [10,11], but still many of their high-temperature properties are either unknown or not systematically explored. One of these is thermal shock behavior.

Thermal shock occurs when a system is subjected to rapid changes in temperature [9]. An abrupt temperature increase gives rise to thermal gradients and hence stress gradients, which may in turn result in damage and catastrophic failure [9]. Therefore, many applications where extreme temperature gradients are required—such as spacecraft propulsion, spacecraft atmospheric entry, immobilization of radioactive waste, combustion, thermoelectric devices, various metallurgical processes, and high-power lasers—are prone to thermal shock [9,12]. Furthermore, physical properties governing the thermal shock behavior are also relevant for other thermomechanical properties such as thermal fatigue [13]. Thermal shock resistance can be described by the thermal shock parameter (R_T), which is defined as

$$R_T=\frac{{\sigma }_f\kappa \left(1-\nu \right)}{\alpha Y},$$

(1)

where σ_f, κ, ν, α, and Y designate flexural strength, thermal conductivity, Poisson’s ratio, the linear coefficient of thermal expansion, and elastic (Young’s) modulus, respectively [9,14]. To increase the thermal shock resistance, R_T should be maximized, which can be achieved by maximizing the numerator and/or minimizing the denominator in Equation (1). Hence, high σ_f is required to enhance resistance to crack propagation, high κ conducts heat away from an active component and minimizes temperature gradients, large ν may give rise to a more ductile response, as well as a combination of low α and small Y leads to a thermal stress reduction [9,14]. It should be noted that an improvement of one physical property in Equation (1) is often accompanied by the deterioration of another property. For instance, a decreased Y is often associated with a larger α value [15], rendering the design of novel thermal shock resistant solids challenging. A typical research strategy is to enhance R_T by concentrating on a single thermomechanical property in Equation (1). For instance, a large R_T value of fused SiO₂ is enabled by very low α [14]. It should also be noted that quantum mechanical predictions, beneficial and efficient in many cases [16], are very challenging for thermal shock behavior since there is an interplay between phonons and electrons in these thermomechanical properties. Hence, to replace the traditional trial-and-error approach—which aims to optimize one of the relevant properties to enhance thermal shock resistance with the knowledge-based design of thermal shock-resistant solids taking into account all (or most) key properties—presents a formidable challenge.

In this work, 14 binary oxides TiO₂, VO₂, CrO₂, MnO₂, NbO₂, RuO₂, RhO₂, TaO₂, OsO₂, IrO₂, SnO₂, PbO₂, SiO₂, and GeO₂ are considered for thermal shock resistance applications using density functional theory [17] in conjunction with the Slack [18] and Debye–Grüneisen model [19]. It should be remarked that many of these thermomechanical properties are affected by extrinsic factors. One of these is microstructure. Flexural strength [20], thermal conductivity [21,22], and the linear coefficient of thermal expansion [23] depend on grain size. For instance, the flexural strength of bulk Ti₃SiC₂, a promising nanolaminate for thermal shock resistance applications, is 330 and 600 MPa at room temperature for grain sizes of 100–200 and 3–5 μm, respectively [20]. The extrinsic factors cannot be straightforwardly addressed by quantum mechanical methods and hence are not considered in the current work. In particular, due to a complex mechanical loading during flexural strength measurements (standard three-point loading) [24], σ_f is not considered in this work. All other four relevant parameters in Equation (1), namely κ, ν, α, and Y, are systematically explored for the common rutile oxides in this study.

2. Methods

Density functional theory [17] was employed in the current work to calculate κ, ν, α, and Y for TiO₂, VO₂, CrO₂, MnO₂, NbO₂, RuO₂, RhO₂, TaO₂, OsO₂, IrO₂, SnO₂, PbO₂, SiO₂, and GeO₂. The Vienna ab initio simulation package (VASP) was used within the framework of the projector augmented wave potentials [25,26,27] and generalized gradient approximation, which were parametrized by Perdew, Burke, and Ernzerhof [28]. The Blöchl correction in the VASP code was applied [29] for these rutile oxides and an integration in the Brillouin zone was carried out by employing the Monkhorst–Pack approach [30] with a k-point mesh of 7 × 7 × 5 (2 × 2 × 2 tetragonal supercell, 48 atoms). The supercell was considered rather than a primitive cell in order to allow for the treatment of diluted alloying (2.1 at.%). No symmetry breaking was observed. Full structural optimization for these tetragonal supercells was made by minimizing the interatomic forces and optimizing the lattice parameters, a and c. The convergence criterion for the total energy was 0.01 meV and a cut-off energy was 500 eV. All configurations were treated as nonmagnetic apart from CrO₂. In the case of CrO₂, spin polarization (ferromagnetic ordering) was also taken into account. Electronic structure analysis of these oxides was performed employing the VESTA software [31].

Two acoustic phonon models were considered. These two phonon models are complementary; one regards the thermal conductivity and the other thermal expansion and temperature-dependent elasticity. Taking Umklapp phonon–phonon scattering into account within the Slack model [18,32], κ values were obtained as

$$\kappa =A\frac{\overline{M}{D}_a^3\delta n^{1/3}}{{\gamma }^{2}T},$$

(2)

where A is a constant, which can be attained as

$$A=\frac{2.43\times {10}^{-6}}{1-0.514{\gamma }^{-1}+0.228{\gamma }^{-2}},$$

(3)

with γ designating the acoustic mode Grüneisen parameter [33,34]. $\overline{M}$ in Equation (2) is the average atomic mass and D_a is the Debye temperature (D) of acoustic phonons calculated as

$$D_a=D{n}^{-1/3}.$$

(4)

Further parameters used in Equation (2), δ³, n, and T are the volume per atom, number of atoms in the 2 × 2 × 2 tetragonal supercell, and absolute temperature, respectively [32]. The values for γ and D_a can be extracted from elastic constants [33,35]. In particular, all tetragonal elastic constants, C₁₁, C₁₂, C₁₃, C₃₃, C₄₄, and C₆₆, were calculated using a method described previously [36]. The tetragonal lattice was strained using a different distortion for each elastic constant (maximum distortion 2%) and the resulting total energy change (quadratic energy–distortion dependence) was utilized to calculate the corresponding elastic constant at 0 K, as detailed in the literature [36]. The elastic constants were also utilized to obtain ν and Y at 0 K within the Hill approximation [37].

The Debye–Grüneisen (acoustic phonon) model was employed for estimating α [19,38,39]. Within the Debye–Grüneisen model [19], the Helmholtz free energy (F) is defined as

$$F=E_{tot}-n{k}_{B}T\left[3{\left(\frac{T}{D}\right)}^3{\int }_0^{D/T}\frac{x^3}{e^x-1}dx-3\ln \left(1-e^{-\frac{D}{T}}\right)-\frac{9D}{8T}\right],$$

(5)

where E_tot and k_B are the total energy at 0 K and the Botzmann constant, respectively. These data were fitted to the Birch–Murnagham equation of state [40] at each temperature to acquire the temperature dependent equilibrium volume and bulk moduli data. The bulk moduli were then used to estimate the temperature dependent Y values, using ν. In the original work on the Debye–Grüneisen model [19], ν was assumed to be constant (0.33) since only metals were considered, but in this work, the ν value was explicitly calculated for each compound minimizing possible errors. In particular, ν was obtained from (3B − 2G)/(6B + 2G)), where B and G are bulk and shear modulus, respectively, within the Hill approximation [37]. The value of α was extracted from the temperature-dependent equilibrium volume datasets. In this work, α was considered instead of volume expansion, since the original Debye–Grüneisen model [19] contains α.

3. Results and Discussion

Figure 1 contains the calculated κ data at 300 K as a function of α for common rutile oxides explored in this study, including TiO₂, VO₂, CrO₂, MnO₂, NbO₂, RuO₂, RhO₂, TaO₂, OsO₂, IrO₂, SnO₂, PbO₂, SiO₂, and GeO₂. Based on the Debye–Grüneisen theory, α and γ are linearly proportional [41]. Hence, an inverse quadratic dependence (parabolic relationship) based on Equation (2) is expected for κ and α. As α increases from 3.3 × 10⁻⁶ to 13.4 × 10⁻⁶ K⁻¹, κ decreases from 35.4 to 1.4 W m⁻¹ K⁻¹. The boundary of the data in Figure 1 is span by SiO₂ (the smallest α value and the largest κ value) and PbO₂ (vice versa). The rest of the data exhibit the values between these two extremes in a parabolic arrangement. The here calculated α and κ value for TiO₂ deviates 29% and 17% from the measured data in the literature [42,43]. It should be mentioned that α is typically within 30% deviated from experiments when the Debye–Grüneisen model is employed [38]. The Slack model is typically less precise, but it reaches a correct order of magnitude [22]. The influence of magnetism is present within these deviations, as probed for CrO₂. For accuracy reasons, the data in Figure 1 are not fitted at this point, but the exact parabolic relationship is discussed below when more datasets are regarded. It should also be noted that the theoretical data obtained by these models may deviate not only due to approximations (single crystals are considered, only acoustic phonons are included in the Slack and Debye–Grüneisen model, no electronic contributions are taken into account within the Slack model, etc.), but also due to difficulties in comparison with available experimental data (polycrystalline samples, impurities, defects, etc.). Nevertheless, important trends are captured in Figure 1. It is also feasible to affect the data in Figure 1 by dilute alloying. Adding 2.1 at.% of Si into TiO₂, increases its κ value by 10.2%. Oppositely, the same amount of Pb leads to a reduction of κ by 15.3%. It should be noted that κ also scales with equilibrium volume in the same fashion as it does with α. This is indicative that the equilibrium volume or bond length governs the thermal response of these isostructural compounds. However, Umklapp phonon–phonon scattering is considered (parabolic relationship between κ and α) as a major physical mechanism for the data shown in Figure 1.

To rationalize the behavior of the boundaries (extremes) in Figure 1, the electronic structure of SiO₂ (the smallest α value and the largest κ value) and PbO₂ (vice versa) is explored in Figure 2. Both Si and Pb are depleted and the majority of charge is attracted by O. This is consistent with ionic bonding. There are important differences between these two rutile phases. While there is essentially no charge localized between Pb and O in PbO₂, a considerable number of electrons are shared by Si and O in SiO₂. Hence, SiO₂ is characterized by ionic–covalent bonding and PbO₂ possesses mainly the ionic contribution to the overall bonding. This is also mirrored in the bond length values for these two extremes, i.e., 1.78 vs. 2.21 Å, respectively. These electronic structure data are consistent with the literature [44,45]. Due to the bonding nature, SiO₂ is expected to have stronger bonds than PbO₂. Stronger bonds thus lead to low α and high κ, as observed in Figure 1 for SiO₂. The opposite occurs for PbO₂.

Figure 3 contains the calculated κ data at 300 K as a function of α for common rutile oxides investigated in this study, previously included in Figure 1, as well as literature values for various semiconductive and insulating phases, including diamond [46], Ge [46], Si [46], S [46], In₂S₃ [47], SiC [48,49], GaN [50,51], Bi₂Te₃ [52], PbTe [13,53], HgTe [54], TiO₂ [42,43], ZnO [55,56], SrTiO₃ [57,58], Kapton (poly-oxydiphenylene-pyromellitimide) [59], and polyvinyl chloride (PVC) [60]. As α increases from 1 × 10⁻⁶ to 61 × 10⁻⁶ K⁻¹, κ decreases from 1000 to 0.1 W m⁻¹ K⁻¹. The κ and α values for the common rutile oxides explored herein are consistent with the functional dependence of the literature data. The obtained inverse square fit for the data in Figure 3 at 300 K gives κ = 763 α⁻², as provided by the solid line, where κ is in units of W m⁻¹ K⁻¹ and α in units of 10⁻⁶ K⁻¹. It is proposed that the constant in the acquired relationship is predominantly determined by the product between the Debye temperature and average atomic mass in Equation (2). These two factors change to a large extent, unlike the other factors, but their product is approximately constant. For instance, this product for isostructural diamond and Ge is 26,443 and 29,322 K u, respectively, where u is the unified atomic mass unit. Furthermore, it is known that small changes in the bonding nature can give rise to a diverse thermal response [61,62]. Since all rutile oxides explored in the current study are isostructural and exhibit similar ionic–covalent bonding (see Figure 2 for extremes), the bond strength is likely the key factor responsible for differences in the thermal properties (Figure 3). For instance, TiO₂ [63] and RuO₂ [64] exhibit an equivalent phonon band structure. With the bond length of 1.97 and 1.99 Å for TiO₂ and RuO₂, respectively, TiO₂ exhibits stronger bonds and should thus possess a higher κ and smaller α value. Indeed, this is observed herein (see Figure 1 for details). Extremes, SiO₂ and PbO₂, undergo the same rationale, as discussed above. To vindicate the whole range of data shown in Figure 3, lattice dynamics of all compounds should be considered. This is beyond the scope of this work, since the Slack and Debye–Grüneisen model is employed instead of full phonon calculations, but the important trends are captured, which is valuable for physical insights and further explorations.

After discussing the behavior of κ and α in Equation (1), the elastic properties are considered. Figure 4 shows the dependence of Y on T for the common rutile oxides explored herein. The Y value slowly decreases with T, as expected, and the difference between the values for these binary oxides is constant. The stiffest rutile oxide is SiO₂ and the softest one is PbO₂, which is in agreement with the bonding analysis (see Figure 2) since SiO₂ possesses the strongest bonds. The obtained elasticity values at 300 K are consistent with the available experimental data, deviating by 7.8% for TiO₂ [65], 2.4% for CrO₂ [2], 16.1% for MnO₂ [66], 2.6% for NbO₂ [67], 5.4% for OsO₂ [68], 9.9% for IrO₂ [69], 0.3% for PbO₂ [70], and 7.9% for SiO₂ [71]. This is acceptable based on the employed exchange-correlation functional, since deviations are commonly within 20% [72]. The calculated internal free parameter for the 4f Wyckoff site (O position) is in a narrow range from 0.345 to 0.348, being 11% deviated from the experiment value [73], but due to the obtained consistency with the elasticity data, this is acceptable. Moreover, the calculated values of Poisson’s ratio are in the range from 0.20 (SiO₂) to 0.32 (PbO₂).

With calculated elasticity (Figure 4) and thermal properties (Figure 1), it is possible to estimate the temperature behavior of R_T, as defined in Equation (1), for the rutile oxides addressed in this study. As argued above, σ_f is obtained under complex mechanical loading and it exhibits the strongest extrinsic response (grain size dependence). Hence, the thermal shock behavior of the common rutile oxides is herein described within density functional theory as R_T/σ_f (units of W m⁻¹ MPa⁻¹ instead of W m⁻¹) and referred to as an intrinsic thermal shock parameter. Figure 5 contains such data as a function of temperature. The intrinsic thermal shock parameter for SiO₂ is the largest in the whole temperature range and decreases from 16.2 to 2.6 W m⁻¹ MPa⁻¹ upon temperature increase from 300 to 900 K. It exceeds the room-temperature value of 11.2 W m⁻¹ MPa⁻¹ measured for nanolaminated Ti₃SiC₂ [74,75] and 7.7 W m⁻¹ MPa⁻¹ for corundum Al₂O₃ [14], which are commonly employed for thermal shock resistance. The rest of the other rutile oxides in Figure 5 exhibit the intrinsic thermal shock parameter lower than 6 W m⁻¹ MPa⁻¹, whereby PbO₂ constitutes the lower bound. Hence, in applications where thermal shock resistant phases are required, Si-based systems are expected to perform well based on the data in Figure 5. At elevated temperatures and under atmospheric conditions, constituting Si is likely to oxidize and enhance the thermal shock resistance. This of course requires an experimental validation, indicating that this work may inspire future investigations. A clue that this notion is correct can be found in the literature on ZrO₂, claiming that additions of SiO₂ enhance its thermal shock resistance [76]. It should also be mentioned that some applications benefit from low κ, such as energy generation and sensing by thermoelectric devices, so that Pb-based systems are of interest, but an interplay with thermal shock resistance and thermal fatigue is challenging to capture and often overlooked in many studies. Holistic approaches, including both primary and secondary properties, are crucial for future design efforts of multifunctional solids.

4. Conclusions

Fourteen rutile oxides have systematically been explored for thermal shock resistance applications using density functional theory in conjunction with the Slack and Debye–Grüneisen model. Four thermomechanical properties (κ, ν, α, and Y) were evaluated, omitting σ_f due to complex mechanical loading and its compelling extrinsic response (grain size dependence). As α increases from 3.3 × 10⁻⁶ to 13.4 × 10⁻⁶ K⁻¹, κ decreases from 35.4 to 1.4 W m⁻¹ K⁻¹, exhibiting a parabolic relationship. The boundary of these data is span by SiO₂ (the smallest α value and the largest κ value) and PbO₂ (vice versa). In a broad comparison with the literature data, an inverse square fit at 300 K was obtained yielding κ = 763 α⁻², where κ is in units of W m⁻¹ K⁻¹ and α in units of 10⁻⁶ K⁻¹. The constant in the acquired relationship may predominantly be determined by the product between the Debye temperature and the average atomic mass. The boundary in the κ–α space for these rutile oxides may be due to stronger bonding in SiO₂ than PbO₂, since SiO₂ is characterized by ionic–covalent bonding and PbO₂ is mainly ionic. This also gives rise to the largest elastic modulus of SiO₂ in the order of 500 GPa at a wide temperature range up to 900 K. Finally, the intrinsic thermal shock resistance is the largest for SiO₂, exceeding some of the competitive phases such as nanolaminated Ti₃SiC₂ and corundum Al₂O₃. It may be argued that at elevated temperatures and under atmospheric conditions, Si-containing systems oxidize, forming SiO₂, and in turn enhance the overall thermal shock resistance.