I. INTRODUCTION
Silicon-based nonoxide structure ceramics, such as Si3N4, SiC, and SiC-matrix composites, have great potential applications in gas turbine engines owing to their superior high-temperature mechanical properties, such as strength and creep resistance. Reference Raj1 However, a major stumbling block in realizing Si-based ceramic turbine hot section components is the recession of these materials in combustion environments because of the volatilization of silica scale by water vapor. Reference Smialek, Robinson, Opila, Fox and Jacobson2,Reference Klemm, Taut and Wötting3 This recession of Si-based ceramics can be prevented by an environmental barrier coating (EBC). Reference Lee4
EBCs are designed to protect Si-based ceramics from corrosion due to water vapor in combustion environments. Therefore, integrity of coatings during the whole duration of service at elevated temperatures is essential. In combustion conditions, flaws caused by collision of foreign particles are unavoidable and the thermal stress arose from mismatched thermal expansion coefficient (TEC) may result in spall and cracking. Reference Dericioglu, Zhu, Kagawa and Kasano5 Thus, the durability of EBCs depends on the tolerance of the coatings to the external mechanical damages, which is closely related to their mechanical properties. Moreover, the thermal stress in the substrate is proportional to Young's modulus and TEC of the given coating material. Reference Lee, Eldridge and Robinson6 To minimize the thermal stress that may lead to undesirable failure and/or creep of the substrate, materials that have the similar TEC with the counterpart substrate and a low Young's modulus are desirable. Reference Lee, Eldridge and Robinson6
Recent experiments found that rare-earth silicates are promising candidates of EBC materials owing to their superior high-temperature stability, long durability in water vapor, and excellent chemical and mechanical compatibility with silicon-based matrix. Reference Lee4,Reference Lee, Fox and Bansal7,Reference Chen, Gao, Liu and Luo8 Among all the studied compounds, ytterbium monosilicate (Yb2SiO5) attracted specific attention since no high-temperature polymorphs have been found for Yb2SiO5. Reference Felsche9 Volatility tests in combustion environments have proved that Yb2SiO5 is nearly inert to the corrosion attack and shows much higher stability than Yb2Si2O7 in the testing conditions. Reference Klemm10 Motivated by its excellent properties, several works focused on the syntheses of Yb2SiO5 have been reported by various approaches, such as solid-state reaction, sol–gel, and hydrothermal syntheses. Reference Chen, Gao, Liu and Luo8,Reference Wen, Dong, He, Wang, Zhou and Zhang11,Reference Khan, Zou, Huang, Fan, Gu, Chen, Zeng, Wang and Cao12 However, reports on the mechanical and thermal properties of Yb2SiO5 are still lacking. Clear understanding of mechanical and thermal properties and correlation of these properties with the structure of Yb2SiO5 can provide valuable insights for the development of rare-earth silicate EBCs.
In this work, mechanical and thermal properties of Yb2SiO5 are investigated from a combination of first-principles calculations and chemical bond theory. The equilibrium crystal structure, second-order elastic constants, and polycrystalline mechanical properties are calculated by first-principles calculations. The anisotropy of elastic behavior is discussed with respect to structural characteristics. Then, the theoretical minimum thermal conductivity and the temperature-dependent thermal conductivity of Yb2SiO5 are predicted from Clarke's model Reference Clarke13,Reference Liu, Wang, Li and Zhou14 and Slack's model, Reference Slack15 respectively. To validate the theoretical calculation, an empirical method based on chemical bond theory is used to calculate bulk modulus and the average linear TEC of Yb2SiO5. Our results demonstrate that Yb2SiO5 is a promising EBC material with a low Young's modulus and low thermal conductivity.
II. COMPUTATION METHODS
Theoretical calculations on structural and elastic properties of Yb2SiO5 were accomplished by the Cambridge Serial Total Energy Package (CASTEP) code, Reference Segall, Lindan, Probert, Pickard, Hasnip, Clark and Payne16 wherein the Vanderbilt-type ultra-soft pseudopotentials Reference Vanderbilt17 were used to represent the interactions between the ionic cores and the valence electrons for Yb, Si, and O atoms. The exchange–correlation energy was treated under the generalized-gradient approximation (GGA). Reference Perdew, Burke and Ernzerhof18 It is known that the additional U correction only strongly affects the electric structures of the rare-earth oxides, whereas the mechanical moduli are insensitive to the localization of partially filled 4f orbitals. Reference Feng, Xiao, Wan, Qu, Huang, Chen, Zhou and Pan19 Thus, the theoretical calculations in this work were mainly based on GGA without consideration of U correction. The plane-wave basis set cutoff energy was fixed at 450 eV in all calculations after convergence test. The special points sampling integration over the Brillouin zone is realized using the Monkhorst–Pack method with special k-points meshes of 2 × 4 × 3. Reference Monkhorst and Pack20
The crystal structure was fully optimized by independently modifying lattice parameters and internal atomic coordinates. The free enthalpy, interatomic forces, and stresses of the unit cell were minimized under the Broyden–Fletcher–Goldfarb–Shanno (BFGS) minimization scheme. Reference Pfrommer, Côté, Louie and Cohen21 The tolerances for geometrical optimization were: differences for total energy within 1 × 10−5 eV/atom, maximum ionic Hellmann–Feynman force within 0.03 eV/Å, maximum ionic displacement within 1 × 10−3 Å, and maximum stress within 0.05 GPa.
The second-order elastic constants were calculated using the strain–stress relationship method implemented by Milman and Warren. Reference Milman and Warren22 To calculate 13 independent elastic constants of Yb2SiO5 (i.e., c 11, c 22, c 33, c 44, c 55, c 66, c 12, c 13, c 23, c 15, c 25, c 35, and c 46), four homogeneous strain patterns were applied and the resulting stresses were calculated by optimizing the internal atomic degrees of freedom. Both positive and negative strains were applied for each strain component, with a maximum strain value of 0.3%. The elastic stiffness was determined from a linear fit of the calculated stress as a function of strain. The criteria for convergence of optimization on atomic internal freedoms were selected as differences in the total energy within 2 × 10−6 eV/atom, Hellmann–Feynman force within 0.006 eV/Å, and maximum ionic displacement within 2 × 10−4 Å. The bulk modulus B and shear modulus G were calculated from the compliance tensor based on Voigt–Reuss–Hill approximation. Reference Voigt23–Reference Hill25 The Young's modulus E H and Poisson's ratio ν were calculated using Hill's bulk modulus B H and shear modulus G H using the following equation: Reference Green26
$$E = {{9{B_{\rm{H}}}{G_{\rm{H}}}} \over {3{B_{\rm{H}}} + {G_{\rm{H}}}}}\quad {\rm{,}}$$
$$\nu = {{3{B_{\rm{H}}} - 2{G_{\rm{H}}}} \over {2\left( {3{B_{\rm{H}}} + {G_{\rm{H}}}} \right)}}\quad {\rm{.}}$$
The minimum thermal conductivity (κmin) of Yb2SiO5 is accessible based on the model proposed by Clarke. Reference Clarke13 He suggested that theoretical minimum conductivity can be calculated after replacing different atoms by an equivalent atom with a mean atomic mass of M/n:
$${\kappa _{\min }} \to {k_{\rm{B}}}{\upsilon _{\rm{m}}}{\left( {{M \over {nd{N_{{\rm{Av}}}}}}} \right)^{ - 2/3}}\quad ,$$
where k
B is the Boltzmann's constant, υm is the average sound velocity, N
Av is the Avogadro's number, d is the density, M is the molecular weight, and n is the number of atoms in the molecule. Clarke used a simplification for υm
$\left( {{\upsilon _{\rm{m}}} = 0.87\sqrt {{E \over d}} } \right)$
without considering anisotropy. Liu et al.
Reference Liu, Wang, Li and Zhou14
proposed a modification to the Clarke's relation by taking anisotropy of sound velocity into account. This modified model has been successfully applied to calculate the κmin of several complex materials, like La2T2O7,
Reference Liu, Wang, Li and Zhou14
MP2O7,
Reference Xiang, Hai and Zhou27,Reference Xiang, Hai and Zhou28
Y4Si2O7N2,
Reference Sun, Liu, Wang, Wang, Zhou and Hu29
and MPO4.
Reference Zhou and Liu30
The average sound velocity υm of polycrystalline material is defined as the following: Reference Anderson31
$${\upsilon _{\rm{m}}} = {\left[ {{1 \over 3}\left( {{2 \over {\upsilon _{\rm{s}}^3}} + {1 \over {\upsilon _{\rm{l}}^3}}} \right)} \right]^{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 3}} \right. \kern-\nulldelimiterspace} 3}}}\quad {\rm{,}}$$
where υl and υs are longitudinal and transverse sound velocities, respectively, which can be determined from the shear modulus G and the bulk modulus B by:
$${\rm{\upsilon }}_{\rm{l}} = \sqrt {{{B + {\raise0.7ex\hbox{${4G}$} \!\mathord{\left/ {\vphantom {{4G} 3}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$3$}}} \over d}} {\rm{,}}\;{\rm{and}}\;{\rm{\upsilon }}_{\rm{s}} = \sqrt {{\raise0.7ex\hbox{$G$} \!\mathord{\left/ {\vphantom {G d}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$d$}}} \quad {\rm{.}}$$
The Debye temperature, ΘD, was calculated according to the following equations: Reference Anderson31
$$\Theta _{\rm{D}} = {h \over {k_{\rm{B}} }}\left[ {{{3n} \over {4\pi }}\left( {{{N_{{\rm{Av}}} d} \over M}} \right)} \right]^{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}} {\rm{\upsilon }}_{\rm{m}} = Cl{\rm{\upsilon }}_{\rm{m}} \quad ,$$
where h is the Plank's constant, k
B is the Boltzmann's constant, n is the number of atoms in the molecular formula, N
Av is the Avogadro's number, d is the theoretical density, and M is the molecular weight.
$l = \left( {{{nd} \over M}} \right)^{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}}$
is defined as the linear density and
$C = {h \over {k_{\rm{B}} }}\left( {{{3N_{{\rm{Av}}} } \over {4\pi }}} \right)^{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}}$
is a constant.
The dependence of thermal conductivity on the temperature is calculated from Slack's model: Reference Slack15
$${\kappa _{\rm{L}}} = A{{\bar M{\Theta ^3}\delta } \over {{\gamma ^2}{n^{{2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-\nulldelimiterspace} 3}}}T}}\quad {\rm{,}}$$
where n is the number of atoms in the primitive unit cell, δ3 is the volume per atom, Θ is the Debye temperature,
$\bar M$
is the average mass of the atoms in the crystal, and A is a physical constant (A ≈ 3.1 × 10−6 if κL is in W m−1 K−1, and δ in Å). The high-temperature limit of the acoustic phonon mode Grüneisen parameter is γ, which is a direct measure of the anharmonicity of phonon, and could be derived from the Poisson's ratio:
Reference Sanditov, Tsydypov and Sanditov32
$$\gamma = {{9\left( {\upsilon _{\rm{l}}^2 - {4 \over 3}\upsilon _{\rm{s}}^2} \right)} \over {2\left( {\upsilon _{\rm{l}}^2 - 2\upsilon _{\rm{s}}^2} \right)}} = {3 \over 2}\left( {{{1 + \nu } \over {2 - 3\nu }}} \right)\quad {\rm{.}}$$
III. RESULTS AND DISCUSSION
A. Structural and chemical bonding properties
The crystal structures of Yb2SiO5 and all its constructive polyhedra are shown in Fig. 1. For Yb2SiO5, there are 64 atoms in one unit cell (Z = 8), which consist of two Yb sites, one Si site, and five O sites. The coordinates of each nonequivalent atomic position are shown in Table I. Four nonequivalent O atoms bond with Si atom to form SiO4 tetrahedra, whereas the “extra” non-silicon-bonded oxygen is surrounded by four Yb atoms in the shape of a slightly distorted tetrahedron. The coordination around Yb atoms is more complicated. Two nonequivalent Yb sites coordinate with six and seven O atoms to form YbO7 and YbO6 polyhedra, respectively. For convenience, we labeled the Yb atom in YbO7 polyhedra as Yb1, and the Yb atom in YbO6 polyhedra as Yb2.
Crystal structure of Yb2SiO5 and typical constructive units in Yb2SiO5.
Theoretical and experimental lattice parameters of Yb2SiO5.
Different electronic exchange–correlation interaction approximation methods [i.e., GGA and localized density approximation (LDA)] are used in the calculation of equilibrium lattice configuration of Yb2SiO5. Compared with the experimental values, Reference Anan’eva, Korovkin, Merkulyaeva, Morozova, Petrov, Savinova, Startsev and Feofilov33 the lattice parameters calculated by GGA are overestimated by 1.63%, whereas those by LDA are underestimated by 2.63%. Therefore, GGA is implemented in the following calculations for its higher accuracy. The calculated lattice parameters of Yb2SiO5 are listed in Table I, together with previously reported experimental data Reference Anan’eva, Korovkin, Merkulyaeva, Morozova, Petrov, Savinova, Startsev and Feofilov33 for comparison. As shown in Table I, the calculated lattice parameters deviate from the experimental data by 0.75% for a, 1.63% for b, and 0.88% for c.
The Mulliken population analysis on Yb–O and Si–O bonds of Yb2SiO5 is conducted, as shown in Fig. 2. In the present work, the population analysis is performed using a projection of the plane-wave states onto a localized basis using a technique described by Sanchez-Portal et al., Reference Sanchez-Portal, Artacho and Soler34 and the resulting projected states are then performed using the formalism described by Segall et al. Reference Segall, Shah, Pickard and Payne35 The bond lengths and Mulliken population of Yb–O and Si–O bonds of another ytterbium silicate, β-Yb2Si2O7, are also included for comparison. Reference Zhou, Zhao, Wang, Sun, Zheng and Wang36 For both silicates, the Mulliken population of Si–O bonds is higher than that of Yb–O bonds, indicating higher level of covalency of Si–O bonds. Moreover, the Mulliken population and bond length of Yb–O bonds are more scattered than those of Si–O bonds, demonstrating severer distortion of Yb–O polyhedra in both materials. For a direct and quantitative view on the distortion of polyhedra, the following equation Reference Carvajal, García-Muñoz, Solé, Gavaldà, Massons, Solans, Díaz and Aguiló37 is implemented to calculate the distortion of polyhedra:
$${\Delta _d} = {1 \over n}{\mathop\sum\limits_n {\left( {{{d{{\left( {{\rm{M}}-{\rm{O}}} \right)}_n} - \left\langle {d\left( {{\rm{M}}-{\rm{O}}} \right)} \right\rangle } \over {\left\langle {d\left( {{\rm{M}}-{\rm{O}}} \right)} \right\rangle }}} \right)} ^2}\quad {\rm{,}}$$
where
$d{\left( {{\rm{M}}-{\rm{O}}} \right)_n}$
and
$\left\langle {d\left( {{\rm{M}}-{\rm{O}}} \right)} \right\rangle$
are the M–O bond length for the nth O atom and the average length of M–O bonds in the polyhedron, respectively. With a higher value of Δ
d
, larger distortion of the polyhedron is expected. The degrees of distortion for all polyhedra constructing Yb2SiO5 and β-Yb2Si2O7 are listed in Table II. The degree of distortion of SiO4 tetrahedron is several orders of magnitude smaller than that of Yb–O polyhedra (i.e., YbO6 and YbO7 units) in both compounds. For Yb2SiO5, the distortion of YO7 unit is larger than that of YO6 unit.
Mulliken population and bond lengths of chemical bonds in ytterbium silicate.
Degrees of distortion of polyhedra of Yb2SiO5 and β-Yb2Si2O7.
Considering the heterogeneity and complexity of chemical bonding in Yb2SiO5, contribution from every species of bonds to the total energy of lattice is different from each other. To quantify this difference, the chemical parameters, such as covalency and bond energy, of each bond are calculated using chemical bond theory introduced in the Appendix. First, based on the analysis of the crystal structure, the chemical formula of Yb2SiO5 is decomposed into the following subformulas:
$$\eqalign{ & {\rm{Yb}}_2 {\rm{SiO}}_5 = {\rm{Yb}}(1){\rm{Yb}}(2){\rm{SiO}}(1){\rm{O}}(2){\rm{O}}(3){\rm{O}}(4){\rm{O}}(5) \cr & \quad \quad \quad \;\; = {\rm{ Yb}}(1)_{1/7} {\rm{O}}(1)_{1/3} + {\rm{Yb}}(1)_{1/7} {\rm{O}}(1)_{1/3}^\prime \quad\quad\quad\quad\quad+ {\rm{Yb}}(1)_{1/7} {\rm{O}}(2)_{1/4} + {\rm{Yb}}(1)_{1/7} {\rm{O}}(2)_{1/4}^\prime \quad\quad\quad\quad\quad+ \;{\rm{Yb}}(1)_{1/7} {\rm{O}}(4)_{1/3} + {\rm{Yb}}(1)_{1/7} {\rm{O}}(5)_{1/4} \quad\quad\quad\quad\quad+ {\rm{Yb}}(1)_{1/7} {\rm{O}}(5)_{1/4}^\prime + {\rm{Yb}}(2)_{1/6} {\rm{O}}(2)_{1/4} \quad\quad\quad\quad\quad+ {\rm{Yb}}(2)_{1/3} {\rm{O}}(3)_{2/3} + {\rm{Yb}}(2)_{1/6} {\rm{O}}(4)_{1/3} \quad\quad\quad\quad\quad+ {\rm{Yb}}(2)_{1/3} {\rm{O}}(5)_{1/2} + {\rm{Si}}_{1/4} {\rm{O}}(1)_{1/3} \quad\quad\quad\quad\quad+{\rm{Si}}_{1/4} {\rm{O}}(2)_{1/4} + {\rm{Si}}_{1/4} {\rm{O}}(3)_{1/3} \quad\quad\quad\quad\quad+ {\rm{Si}}_{1/4} {\rm{O}}(4)_{1/3} \quad {\rm{,}} \cr}$$
the superscript in the formula marks the oxygen atoms having the same crystallography site with different bond lengths in one polyhedron. Then, the desired parameters of the subformula are calculated by Eqs. (A1)-(A8) and are listed in Table III. It is worth mentioning that the bond length for each bond is from the experimental data Reference Anan’eva, Korovkin, Merkulyaeva, Morozova, Petrov, Savinova, Startsev and Feofilov33 rather than the theoretical results from density functional theory calculations. The calculated covalency (f c μ) of Si–O bonds is much higher than that of Yb–O bonds, which is consistent with the Mulliken analysis on the bonds. In addition, the bond energies of Si–O bonds are also much higher than that of Yb–O bonds, showing strong bonding nature of Si–O bonds when compared with Yb–O bonds. Meanwhile, the bond energies of Yb–O bonds in YbO6 unit are considerably higher than those in YbO7 unit, indicating that adding an “extra” O atom to the octahedral structure decreases the bonding strength of original Yb–O bonds efficiently. This can also be verified simply by Pauling's rules Reference Pauling38 : the strength of a bond equals the ionic valence divided by the coordination number, which yields a strength of 0.43 for Yb–O bonds in YbO7 polyhedra, 0.5 for Yb–O bonds in YbO6 octahedra, and 1 for Si–O bonds.
Chemical parameters of subformulas for Yb2SiO5 by chemical bond theory.
In general, Yb2SiO5 shows a high level of distortion in the structure, and the bonding nature is quite complex and heterogeneous. SiO4 tetrahedra are strongly constructed and show rigid characteristic, whereas Yb–O units are relatively weak and soft. These structural features play an essential role in the elastic and thermal properties of this material, which is discussed in the following sections.
B. Elastic properties
The calculated elastic constants are listed in Table IV, accompanied by those of β-Yb2Si2O7 as comparison. To the authors' knowledge, the full sets of elastic constants of ytterbium silicate are not available from experiments owing to the low symmetry of monoclinic structure. The mechanical stability of a crystal structure can be justified from the Born–Huang criterion. Reference Born and Huang39 For a monoclinic structure, the conditions are quite complex and have been omitted here for simplicity (one may refer to Ref. Reference Wu, Zhao, Xiang, Hao, Liu and Meng40 for full expressions of the criterion). The criteria are completely satisfied by Yb2SiO5, demonstrating its stability under elastic strain perturbations. Among all constants, c 11, c 22, and c 33 represent the stiffness against principal strains, whereas c 44, c 55, and c 66 correspond to resistance to shear deformations. A common feature for both materials is that c 11, c 22, and c 33 are much higher than c 44, c 55, and c 66. Moreover, these six constants of Yb2SiO5 are obviously lower than those of β-Yb2Si2O7, indicating less resistance of Yb2SiO5 to the principal and shear deformations.
Calculated independent second-order elastic constants c ij (in GPa) for Yb2SiO5 and β-Yb2Si2O7 ceramics.
The calculated mechanical properties of Yb2SiO5 are calculated and tabulated in Table V. All the elastic parameters of Yb2SiO5 are lower than those of β-Yb2Si2O7. The bulk moduli of these two compounds are rather small, showing that they are not relatively incompressible. Unfortunately, no experimental results are available. Therefore, to verify our first-principles calculations, the bulk modulus of Yb2SiO5 is also estimated using chemical bond theory introduced in the Appendix. The results are listed in Table VI. It is seen that the bulk moduli of Si–O bonds are larger than those of Yb–O bonds, showing less compressibility of Si–O bonds. The bulk modulus for Yb2SiO5 polycrystalline is 110 GPa, which agrees with the value obtained by first-principles calculations quite well. Consistence of bulk moduli from different approaches ensures the reliability of our calculations.
Calculated bulk modulus B H (GPa), shear modulus G H (GPa), Young's modulus E (GPa), Poisson's ratio (ν), and Pugh's ratio (G/B) of Yb2SiO5 and β-Yb2Si2O7 ceramics.
Estimation of bond volume νμ (Å3), lattice energy density u μ (GPa), bulk modulus B (GPa), and linear TEC α (×10−6 K−1) of Yb2SiO5 from chemical bond theory.
The shear modulus, G, describes the resistance of a material against a shape change. Interestingly, for both materials, the shear moduli are substantially smaller than bulk moduli, as shown in Table V. The G/B ratio (namely Pugh's ratio) is 0.542 and 0.448 for Yb2SiO5 and β-Yb2Si2O7, respectively. It is widely accepted that the ductility of a ceramic is estimated by the ability of performing shear deformation, which is quantified by Pugh's ratio. For a ductile material, the Pugh's ratio is no larger than 0.541. The Pugh's ratio for brittle ceramics is 0.85 for SiC Reference Lambrecht, Segall, Methfessel and Schilfgaarde41 and 0.76 for ZrC Reference Wu, Chen, Struzhkin and Cohen42 ; whereas for the well-known damage-tolerant layered ternary carbides (the so-called MAX phases), they are less than 0.5. Reference Wang and Zhou43 Small values of G/B suggest that Yb2SiO5 and β-Yb2Si2O7 possess low shear deformation resistance and good machinability, which result in damage tolerance in experiments. The damage tolerance of β-Yb2Si2O7 has been confirmed by Zhou et al., Reference Zhou, Zhao, Wang, Sun, Zheng and Wang36 whereas the experimental investigation on Yb2SiO5 is still in dearth.
The intrinsic hardness of a material is a highly complex property, which is difficult to be described by a formal theoretical definition. Chen et al. Reference Chen, Niu, Li and Li44 proposed an empirical model, which not only correlates with shear modulus G but also with bulk modulus B. In their model, the Vickers hardness is estimated by the following formula:
$${H_{\rm{V}}} = 2{\left( {{k^2}G} \right)^{0.585}} - 3\quad {\rm{,}}$$
where k is the Pugh's modulus ratio. This model is performed well across a large class of materials and hardness values. Reference Chen, Niu, Li and Li44 Here, we use this model to estimate the hardness of Yb2SiO5, which yields only 6.9 GPa. Relatively soft characteristics of Yb2SiO5 result from a small shear modulus, which originates from the weak bonding between Yb and O atoms, as in the case of MAX phases. Reference Wang and Zhou43
The Young's modulus, E, measures the resistance of materials against uniaxial tension. The calculated Young's modulus of polycrystalline Yb2SiO5 is 132 GPa, which is 21% smaller than that of β-Yb2Si2O7. As shown in Table V, the Young's moduli along three principal directions are different, showing the anisotropy in elastic behavior of Yb2SiO5. To obtain a clear and complete representation of the elastic anisotropy of Yb2SiO5, the variation of Young's modulus as a function of crystal orientation is necessary. The direction-dependent Young's modulus, E, for Yb2SiO5 is given by the following equation: Reference Nye45
$$\eqalignb{ & {1 \over E} = l_1^{\,4}{s_{11}} + 2l_1^2l_2^2{s_{12}} + 2l_1^2l_3^2{s_{13}} + 2l_1^3{l_3}{s_{15}} + l_2^{\,4}{s_{22}} \quad\quad+ 2l_2^2l_3^2{s_{23}} + 2{l_1}l_2^2{l_3}{s_{25}} + l_3^{\,4}{s_{33}} + 2{l_1}l_3^3{s_{35}} \quad\quad+ l_2^2l_3^2{s_{44}} + 2{l_1}l_2^2{l_3}{s_{46}} + l_1^2l_3^2{s_{55}} + l_1^2l_2^2{s_{66}}\quad {\rm{,}} \cr}$$
where s ij is the elastic compliance and l 1, l 2, and l 3 are the directional cosines of angles with the three principal directions, respectively.
The surface contour of the Young's modulus of Yb2SiO5 is shown in Fig. 3(a), and the planar projections of Young's modulus for (100), (010), and (001) crystallographic planes are also shown in Fig. 3(b). For different crystallographic planes, A and B directions represent different crystallographic directions: A direction represents [100] and B represents [010] for (001) plane. In the case of (010) plane, they are [100] and [001] directions; whereas for (100) plane, they are [001] and [010] directions. The elastic anisotropy of Yb2SiO5 is clearly illustrated. The anisotropy of Young's modulus on (010) plane is stronger than that on the other two planes. The maximum and minimum values of Young's modulus are deflected from the principal axes. As suggested by Mohapatra and Eckhardt, Reference Mohapatra and Eckhardt46 the anisotropy in elastic moduli is mainly affected by the off-diagonal elements of the compliance matrix. The degree of anisotropy of elastic moduli decreases significantly if the off-diagonal elements (i.e., s 12, s 13, s 23, s 15, and s 35 in this case) are omitted in the calculation of Young's moduli along different crystallographic directions. For simplicity, the procedure is not presented here. Using Eq. (12), the maximum and minimum Young's moduli are 173 and 102 GPa. The high- and low-modulus directions correspond to the high- and low-fracture energy directions. Figure 3 shows the anisotropy of the elastic properties of Yb2SiO5 and also the expected low- and high-fracture energy direction in this material. With this information, the most important directions for the mechanical property measurements and applications are clarified.
(a) Surface contour of direction-dependent Young's modulus of Yb2SiO5 and (b) planar projections on (100), (010), and (001) crystallographic planes.
C. Thermal properties
TEC, resulting from the anharmonic vibration of the lattice at a finite temperature, is a fundamental parameter. Theoretical evaluation of TEC is a challenge for material scientists. Several models have been proposed to estimate TEC: quasi-harmonic approximation, Reference Blanco, Francisco and Luaña47 molecular dynamics (MD), Reference Car and Parrinello48 chemical bond theory, Reference Zhang, Li, Zhou and Pan49 and so on. The chemical bond theory proposed by Zhang et al. Reference Zhang, Li, Zhou and Pan49 is an effective method to predict the linear expansion coefficients of complex crystals based only on the crystal structure data. The calculated results show excellent agreement with the experimental values. Reference Zhang, Li, Zhou and Pan49 Therefore, the linear expansion coefficient of Yb2SiO5 is evaluated by chemical bond theory in our work. The calculated results are listed in Table VI. Considering that the bond strength of Yb–O bonds is weaker than that of Si–O bonds, it is not surprising that the contributions from Yb–O bonds are quite larger than those from Si–O bonds. The expansion coefficient of Yb2SiO5 is 5.24 × 10−6 K−1, theoretically. Experimentally, however, TEC of Yb2SiO5, which is quite rare and scattered, extends from 3.5 × 10−6 to 7.1 × 10−6 K−1. Reference Lee, Fox and Bansal7,Reference Klemm10,Reference Khan, Zou, Huang, Fan, Gu, Chen, Zeng, Wang and Cao12 The discrepancy between experimental and theoretical data shows the lack of knowledge on this material and vividly illustrates that more delicate investigations are needed on TEC of Yb2SiO5, experimentally and theoretically. To accurately calculate the TEC of a material, one may refer to ab initio MD, which will be considered in our future works.
The behavior of thermal transportation at different temperature ranges is an important factor that needs to be considered when selecting an applicable EBC. As shown in Eqs. (3) and (7), estimation of the intrinsic thermal conductivity of Yb2SiO5 relies on the knowledge of sound velocities (υl, υs, and υm), Debye temperature (ΘD), and Grüneisen parameter (γ). The sound velocities and Debye temperature are derived from obtained elastic moduli and density of equilibrium structure by Eqs. (4)–(6). The calculated parameters for Yb2SiO5 are listed in Table VII, together with those of β-Yb2Si2O7 for comparison. It is reasonable that the sound velocities of Yb2SiO5 are significantly lower than those of β-Yb2Si2O7 because smaller values of elastic moduli and higher value of density of Yb2SiO5 result in smaller sound velocities. Relatively low Debye temperature of Yb2SiO5 also suggests relatively weak atomic vibrations in the lattice compared with β-Yb2Si2O7.
Theoretical density d (in g cm−3), sound velocities (υl, υs, υm in km s−1), Debye temperature ΘD (in K), Grüneisen parameter γ, and minimum thermal conductivity (in W m−1 K−1) of polycrystalline Yb–Si–O compounds.
The intrinsic lattice thermal conductivity κL of a material is determined by different phonon scattering mechanism in different temperature ranges. Reference Clarke13 At a relatively high temperature, the dominant mechanism of phonon scattering is the Umklapp processes, in which the acoustic phonon branches interact with each other to transport heat. The Slack's model [i.e., Eq. (7)] Reference Slack15 is a suitable approach to describe the temperature-dependent thermal conductivity when the Umklapp scattering is dominant. The temperature-dependent thermal conductivity of Yb2SiO5 estimated from the Slack's model is shown in Fig. 4. With increase in the temperature, the thermal conductivity of Yb2SiO5 decreases over the studied temperature range. Then, if the temperature increases further, the phonon mean-free path decreases to the average atomic distance, and thus the thermal conductivity will approach its minimum. Reference Clarke13 The minimum thermal conductivity κmin is evaluated by the modified Clarke's model, as illustrated by Eqs. (3) and (4). The minimum thermal conductivity of Yb2SiO5 is predicted to be 0.74 W m−1 K−1. The theoretical value is slightly lower than the measured one (1.3–1.4 W m−1 K−1). Reference Jacobson, Fox, Smialek, Opila, Dellacorte, Lee, Cramer and Covino50 The discrepancy may either derive from the experimental uncertainty or the infrared radiation at high temperatures. Combining these two models, the behavior of thermal transportation of Yb2SiO5 can be described as: with the increase of temperature, the thermal conductivity of Yb2SiO5 declines as κ = 893.82/T; while, at around 1200 K, it reaches a minimum value and remains unchanged with temperature.
Temperature dependence of thermal conductivity of Yb2SiO5. The minimum thermal conductivity (red dash line) is also shown.
The theoretical minimum thermal conductivities are 0.74 and 1.12 W m−1 K−1 for Yb2SiO5 and β-Yb2Si2O7, respectively. Extremely low thermal conductivity could be traced back to the characteristic of the structure. Zhou et al. Reference Zhou, Zhao, Wang, Sun, Zheng and Wang36 concluded that the low thermal conductivity of β-Yb2Si2O7 is mainly attributed to the heterogeneous bond strength. Rigid SiO4 units are efficient in phonon transport, whereas weakly bonded YbO6 units intensify the scattering of phonon. These structure features are also applicable to Yb2SiO5. Moreover, as discussed in Sec. III. A, adding an “extra” O atom to the YbO6 unit decreases the bonding strength and ordering of original octahedron, which would aggravate the phonon scattering, decrease the phonon transportation, and result in much lower thermal conductivity of Yb2SiO5.
For a desirable EBC material, the stability in the hot gas atmosphere is essential. Meanwhile, certain mechanical properties are also required to maintain the integrity of coating and minimize the thermal stress. The theoretical investigations on the mechanical and thermal properties of Yb2SiO5 in this work ascertain that it is a promising thermal and environmental barrier coating material (T/EBC). Small G/B indicates intrinsic quasi-ductility and good damage tolerance, which ensures the integrity of the coating under the impact of foreign particles and thermal cycles and avoids the cracking derived from the mismatched TEC. Compared with other silicates (e.g., 168 GPa for β-Yb2Si2O7, Reference Zhou, Zhao, Wang, Sun, Zheng and Wang36 155 GPa for γ-Y2Si2O7, Reference Sun, Zhou, Wang and Li51 and 150 GPa for Y2SiO5 Reference Luo, Wang, Wang, Li and Hu52 ), the Young's modulus of Yb2SiO5 is fairly small, which benefits in minimizing the thermal stress on the substrate, and increases the thermal shock resistance of Yb2SiO5. Moreover, the theoretical minimum thermal conductivity of Yb2SiO5 is substantially small, which endues it a promising candidate as T/EBC.
IV. CONCLUSIONS
Bonding characteristics, elastic stiffness, and thermal conductivity of Yb2SiO5 are investigated by first-principles total energy calculations and chemical bond theory. The theoretical results reveal the heterogeneous bonding nature of Yb2SiO5, i.e., Si–O bonds are stronger than Yb–O bonds. The second-order elastic constants and mechanical properties of Yb2SiO5 are predicted. The low G/B ratio suggests that Yb2SiO5 is a quasi-ductile ceramic. The elastic anisotropy of Yb2SiO5 is discussed by calculating the direction-dependent Young's modulus of Yb2SiO5. The Debye temperature of Yb2SiO5 is determined to be 380 K. Then, the elastic moduli, sound velocities, and Debye temperature are used to calculate the temperature-dependent thermal conductivity of Yb2SiO5. The intrinsic thermal conductivity of Yb2SiO5 decreases as 893.82/T and approaches a minimum value (0.74 W m−1 K−1) at high temperatures. Our theoretical results show that weak Yb–O bonds play an important role in the mechanical and thermal properties of Yb2SiO5 and highlight the potential of Yb2SiO5 as a promising candidate for thermal and environmental barrier coating.
ACKNOWLEDGMENTS
The authors are grateful to Dr. Bin Liu from Oak Ridge National Laboratory, USA, for the English revision. This work was supported by the National Outstanding Young Scientist Foundation for Y.C. Zhou under Grant No. 59925208 and the Natural Sciences Foundation of China under Grant No. 50832008.
APPENDIX: Chemical Bond Theory
Chemical bond theory, developed by Phillips and Van Vechten, Reference Phillips and Van Vechten53 Van Vechten, Reference Van Vechten54 Levine, Reference Levine55 and Xue and Zhang, Reference Xue and Zhang56 provides an efficient method to study the properties of the chemical bonds in complex crystals and to predict the physical properties of a crystal from a bonding standpoint. The implementation of the theory relies on the decomposition of complex crystal into a linear combination of subformula of various binary crystals. Reference Xue and Zhang56
According to the chemical bond theory, for a giving chemical formula of a complex crystal, it is always possible to decompose it into the following simple subformula: Reference Xue and Zhang56
$$A_{a1}^1A_{a2}^2 \cdots A_{ai}^iB_{b1}^1B_{b2}^2 \cdots B_{bj}^{\,j} = \mathop\sum\limits_{i,j} {A_{mi}^iB_{nj}^{\,j}} = \mathop\sum\limits_\mu {{{\left( {{A_m}{B_n}} \right)}^\mu }} \quad ,$$
$$mi = {{N\left( {{B^{\,j}} - {A^i}} \right) \times ai} \over {{N_{\rm{C}}}\left( {{A^i}} \right)}},\;nj = {{N\left( {{A^i} - {B^{\,j}}} \right) \times bj} \over {{N_{\rm{C}}}\left( {{B^{\,j}}} \right)}}\quad {\rm{,}}$$
where N C(A i ) and N C(B j ) are the coordination numbers of A i and B j ions in the crystal. N(B j − A i ) is the nearest coordination fraction contributed by A i ion and vice versa. (A m B n )μ is one of subformulas. Since the fictitious binary crystal contains only one species of chemical bond, its properties can be easily estimated.
For a fictitious binary crystal (A m B n )μ, the total lattice energy U μ can be separated into the ionic part U i μ and the covalent part U c μ, which can, respectively, be expressed as: Reference Liu, Zhang and Wu57
$$U_{\rm{i}}^\mu = {{1270\left( {m + n} \right){Z^ + }{Z^ - }} \over {{l^\mu }}}\left( {1 - {{0.4} \over {{l^\mu }}}} \right)f_{\rm{i}}^\mu \;\left( {{\rm{kJ}}\;{\rm{mo}}{{\rm{l}}^{ - 1}}} \right)\quad ,$$
$$U_{\rm{c}}^\mu = 2100m{{{{\left( {{Z^ + }} \right)}^{1.64}}} \over {{{\left( {{l^\mu }} \right)}^{0.75}}}}f_{\rm{c}}^\mu \;\left( {{\rm{kJ}}\;{\rm{mo}}{{\rm{l}}^{ - 1}}} \right)\quad {\rm{,}}$$
where Z + and Z − are the valence of cation and anion in the binary crystal and l μ is the length of μ type bond. f i μ and f c μ are the ionicity and covalency of A–B bond, and are defined as: Reference Levine55
$$f_{\rm{i}}^\mu = {{{{\left( {{C^\mu }} \right)}^2}} \over {{{\left( {E_{\rm{g}}^\mu } \right)}^2}}},f_{\rm{c}}^\mu = {{{{\left( {E_{\rm{h}}^\mu } \right)}^2}} \over {{{\left( {E_{\rm{g}}^\mu } \right)}^2}}}\quad {\rm{,}}$$
$$f_{\rm{i}}^\mu + f_{\rm{c}}^\mu = 1\quad {\rm{,}}$$
where E g μ is the average energy bond gap and is composed of homopolar Eh μ and heteropolar C μ parts. These two parts of energy can be estimated by: Reference Levine55
$$E_{\rm{h}}^\mu = {{39.74} \over {{{\left( {{l^\mu }} \right)}^{2.48}}}}\;{\rm{(eV)}}\quad {\rm{,}}$$
$${C^\mu } = 14.4{b^\mu }\;{\rm{exp}}\left( { - k_{\rm{s}}^\mu {r^\mu }} \right)\left( {{{{{\left( {Z_{\rm{A}}^\mu } \right)}^ * }} \over {{r^\mu }}} - \left( {{n \mathord{\left/ {\vphantom {n m}} \right. \kern-\nulldelimiterspace} m}} \right){{{{\left( {Z_{\rm{B}}^\mu } \right)}^ * }} \over {{r^\mu }}}} \right)\;\quad\quad({\rm{eV}}){\rm{ }}({\rm{if}}\;n > m)\quad {\rm{,}}$$
$${C^\mu } = 14.4{b^\mu }\;{\rm{exp}}\left( { - k_{\rm{s}}^\mu {r^\mu }} \right)\left( {\left( {{m \mathord{\left/ {\vphantom {m n}} \right. \kern-\nulldelimiterspace} n}} \right){{{{\left( {Z_{\rm{A}}^\mu } \right)}^ * }} \over {{r^\mu }}} - {{{{\left( {Z_{\rm{B}}^\mu } \right)}^ * }} \over {{r^\mu }}}} \right)\;\quad\quad({\rm{eV}})\;({\rm{if }}m > n)\quad {\rm{,}}$$
where b μ is a structural correction factor, r μ = l μ/2 is the average ion radius expressed in angstroms, exp(−k s μ r μ) is the Thomas–Fermi screening factor, and (Z A μ)* and (Z B μ)* are the effective valence electron numbers of A and B ions, respectively.
The lattice energy density u μ of a binary crystal, which is an important parameter in estimating the bulk modulus and linear TEC of a binary crystal, is defined as: Reference Zhang, Li, Li, Zhou and Cao58
$${u^\mu } = {{{U^\mu }} \over {{N_{{\rm{Av}}}}{n^\mu }v_{\rm{b}}^\mu }}\quad ,$$
where N Av is the Avogadro constant, n μ is the number of the chemical bond in one formula unit, and v b μ is the volume of μ type chemical bond. Then, the bulk modulus B μ of a certain type of binary crystals can be calculated on the basis of lattice energy density: Reference Zhang, Li, Li, Zhou and Cao58
$${B^\mu } = {\delta ^\mu } + {{{u^\mu }} \over {{\beta ^\mu }}}\quad ,$$
where δμ is a constant and βμ is a proportion factor depending on the average valence and average coordination number of subformula A m B n . The bulk modulus of the complex crystals is expressed as: Reference Zhang, Li, Li and Zhou59
$${B_{\rm{m}}} = {1 \over {{k_{\rm{m}}}}}\quad ,$$
$${k_{\rm{m}}} = {{\Delta V} \over {{V_{\rm{m}}}\Delta P}} = {1 \over {{V_{\rm{m}}}}}\mathop\sum\limits_\mu {{{\Delta {V_\mu }} \over {\Delta {P_\mu }}}} = {1 \over {{V_{\rm{m}}}}}\mathop\sum\limits_\mu {{V^\mu }{k^\mu }} = {1 \over {{V_{\rm{m}}}}}\mathop\sum\limits_\mu {{{{V^\mu }} \over {{B^\mu }}}} \quad ,$$
where B μ is the bulk modulus of μ type chemical bond, and V μ is the volume of μ type bond in one molecule.
The linear TEC of a complex crystal is closely related to the lattice energy and can be estimated by the following equations: Reference Zhang, Li, Zhou and Pan49
$$\alpha = \mathop\sum\limits_\mu {{F^\mu }{\alpha ^\mu }} \quad ,$$
$${\alpha ^\mu } = - 3.1685 + 0.8376{\chi ^\mu }\quad ,$$
$${\chi ^\mu } = {{{k_{\rm{B}}}Z_{\rm{A}}^\mu N_{{\rm{CA}}}^\mu } \over {{U^\mu }{\Delta ^\mu }}}{\beta ^\mu },\;{\beta ^\mu } = {{m\left( {m + n} \right)} \over {2n}}\quad ,$$
where parameter Δμ is a correction parameter from the analysis of experimental results, which depends on the position of cation in the periodic table. F μ is the fraction of μ type bond in the complex crystal.
