Characterizing Property of States: Effect of Defects on the Coefficient of Thermal Expansion and the Specific Heat Capacity of ZrB₂

Jude O. Ighere, P. Alex Greaney

School of Mechanical, Industrial & Manufacturing Engineering, Oregon State University, Corvallis, OR, USA.

DOI: 10.4236/njgc.2020.102002                

Abstract

Thermal storage potential and thermal expansion are characteristic properties for extreme applications. ZrB₂ is a candidate for advanced applications in aircraft and fusion reactors. This article presents density functional theory calculations of its states, microstructure and quasi-harmonic levels calculations of thermophysical properties. Band structure highlighted dynamical instability with metallic impurities in ZrB₂ structure based on frequency modes. The observed projected density of states (PDOS) appropriate 4d orbital of Zr dominated at low frequency both in perfect crystal in the presence or absence of covalent impurities while B 2s and 2p orbitals dominate higher frequency states. Temperature dependency and anisotropy of coefficient of thermal expansion (CTE) were evaluated with various impurities. Various thermodynamic properties like entropy and free energy were explored for degrees of freedom resulting from internal energy changes in the material. Computed results for heat capacity and CTE were compared to available numerical and experimental data.

Keywords

Zirconium Diboride, ZrB₂, Thermal Expansion, Defects, Molecular Dynamics, Specific Heat Capacity

1. Introduction

Zirconium diboride (ZrB₂) attracts much attention as it is used in various high temperature applications including nuclear reactors, turbine engines and leading-edge aircraft due to its unique thermal properties at extreme conditions. This ultra-high temperature ceramics (UHTC) reportedly melts at the range of 3000˚C - 3245˚C [1] [2]. The characteristic properties of expansion and heat capacity are therefore critical to inform life design methods for thermal stress optimization. This article is concerned with the impact of defects on its thermal properties and the hexagonal lattice structure of ZrB₂. With relatively low density, ZrB₂ is a list of few candidates for hypersonic flight which generates high temperatures at its leading edges.

Its atoms do not deviate from their lattice sites hence the thermal and mechanical stability which is ascribed to its well layered in-plane and out-of-plane bonding in the hexagonal honeycomb (Figure 1). The behavior of the material under thermal stress is quantified using heat capacity and thermal expansion. However, the perfect crystal of ZrB₂ rarely exists due to its partially covalent nature, and induced defects from synthesis and sintering processes [3] [4] [5]. Defects such as surface oxides, metallic impurities and carbides are unavoidably created during processing stages or exist in nature. These impurities, when present, result in microscopic structural changes, and hence redistribution in localized electron density. This induced variation in the mechanical reliability of the material creates phonon or electron (carrier) scattering differences evident in its thermal properties or processes.

First principle calculations are critical at such microscopic levels to contrast the perfect ZrB₂ crystals to one with impurities and correlate electrical and thermodynamic properties. In this article, band structures and density of states (DOS) were investigated and the defect commonalities to heat capacity and thermal expansion explored.

1.1. Coefficient of Thermal Expansion of ZrB₂

In considering ZrB₂ as the leading favorite in the internal design of nuclear reactors, one should critically evaluate its lifecycle potential based on expansivity [6] [7]. Studies have shown that the thermal expansion of ZrB₂ is anisotropic along the various lattice directions [8]. This directionally dependent property is attributed to the varying bond strengths along the different axes which determine stiffness and therefore the expansion [9] [10]. Researchers have attributed

Figure 1. Hexagonal honeycomb layers of boron-boron atoms alternate with layers of Zirconium atoms centered in the hexagons.

the characteristic stiffness of ZrB₂ crystal structure to the strong in-plane boron-boron sigma covalent bonds (Figure 1). [9] [11]. Generality of the expression for coefficient of expansion holds true on the assumptions: (i) that the expansion is insensitive to the individual axial-direction (ii) the individual molecules do not expand but merely vibrates due to kinetic energy. Therefore, for ZrB₂, the essence of linear expansion is also considered to quantify the anisotropic thermal expansion especially in the z-direction. The coefficient of linear thermal expansion is given as ${\alpha }_{x,y,z}=\left(1/L_i\right)\left(\text{d}L/\text{d}T\right)$, where L_i represents specific axial-length (x, y, z) at temperatures T and initial.

For a typical solid, the coefficient of thermal expansion (CTE) (or volumetric expansivity) is generally expressed as:

$\alpha =\frac{V_f-V_i}{V_i\left(T_f-T_i\right)}$ (1)

where $\alpha$ is the volumetric expansion coefficient, $V_{f,i},T_{f,i}$ are the respective initial and final volumes and temperatures. However, in potentially anisotropic materials, the linear (directional) expansion along the different directions can be calculated using:

$\alpha \left(l\right)=\frac{1}{l}\left(\frac{\delta l\left(T\right)}{\delta T}\right)$ (2)

where l dimension of the crystal lattice in any direction x, y, z and $\alpha \left(l\right)$ is respective linear thermal expansion.

In leading edge applications, the sharp-noses (Figure 2) of supersonic aircraft moving through space with speed many times the speed of sound are subjected to extremely high thermal stress. Uneven expansion and contraction across the temperature gradient poses a threat to structural reliability of the material. Thus, this work investigated temperature dependency of CTE for ZrB₂ and the impact of impurities.

1.2. Specific Heat Capacity of ZrB₂

Several studies have focused on thermal conductivity and diffusivity of ZrB2-based ceramics but only a few have given attention specific heat capacity.

Figure 2. Depicting temperature gradient and regional component in ramping temperatures [12].

This property is essential in characterizing the state of the microstructures of materials for improved tailoring [13] [14]. It can be defined at constant volume or pressure and correlated respectively to the internal energy or enthalpy of the crystal structure.

Phonon-based (lattice) constant volume heat capacity was formulated by Einstein and later Debye as shown in Equations (3) and (4) respectively:

$C_{Ev}=3N{K}_B\frac{x^2{\text{e}}^x}{{\left({\text{e}}^x-1\right)}^2}$ (3)

$C_{Dv}=9N{K}_{B}T{\left(\frac{T}{{\theta }_D}\right)}^3{\displaystyle {\int }_0^{{\theta }_D/T}\left(\frac{x^4{\text{e}}^x}{{\left({\text{e}}^x-1\right)}^2}\right)\text{d}x}$ (4)

where ${\theta }_D$ and N represent the Debye temperature and the Avogadro and $x=\hslash \omega /k_{B}T$. This relationship is simplified at low and can be written as:

$C_v=\beta T^3$ (5)

The coefficient, $\beta$ (specific heat coefficient) is correlated to the Debye temperature ${\theta }_D\left(0\right)$ via:

${\theta }_D\left(0\right)={\left(\frac{1944}{\beta }\right)}^{\frac{1}{3}}$ (6)

Both models show that the transfer of phonons in lattices contributes to measurable change in the energy. And Debye’s model further attributes the interaction between harmonic oscillating atoms to the exponential increase in heat capacity with temperature. Also, analyzed in this article was analyzed the free energy, entropy and enthalpy of the structures in a bid to explain the heat capacity and internal energy.

2. Computational Details

The density function theory (DFT) calculations for this study were performed using Vienna Ab-initio Simulation Package (VASP) [15] [16]. The projected augmented wave (PAW) exchange potential was used in combination with the correlation functional developed by Perdew, Burke and Ernzerhof (PBE). A cutoff of 600 eV was used for the kinetic energy for the planewave, and a Monkhorst-Pack k-point mesh of 4 × 4 × 4. The free (total) energy change and band structure energy change convergence point was set to less than 1 × 10⁻⁴ eV between steps. For a more accurate band structure and density of states, higher k-point grid 12 × 12 × 12 was further used. A classical molecular dynamic simulation data was also employed for comparison. Phonopy (open source package) was used for phonon calculations. Package uses statistical thermodynamic expressions to compute free energy (F), heat capacity (C), entropy (S), and enthalpy (H). Smearing method was used via Phonopy to calculate the phonon density of states (DOS) on a sampling mesh.

The classical MD simulations were performed using Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) to calculate the coefficient of thermal expansion (CTE) of perfect ZrB₂ and contributions from impurities [17]. Tersoff interatomic potential was used to describe the atomic interactions [18]. Tersoff potential has been used previously to model the properties of Zr, B, Hf, ZrB₂ and SiC [5] [10] [18]. Periodic boundary conditions were applied to model the ZrB₂ crystal structure to enable interaction across boundary. Langevin thermostat was used to control the temperature within NPT ensemble. A timestep of 5 fs was employed with 0.6 ns equilibration time for within a simulation duration of 15 ns. Literature valures for ZrB₂ lattice constant were adopted, $a_o,b_o=3.035$ Å and $c_o=3.223$ Å. All atoms with same z-coordinates are termed same plane, with alternating B and Zr layer. Added validation procedure to ensure that the crystal lattice generated (see Figure 3(b)) does not deviate was applied and VESTA was used to load and generate XRD patterns. Figure 3 compares calculated to published experimental XRD patterns for perfect ZrB₂ [8].

Figure 3. (a) XRD pattern for perfect ZrB₂ crystal lattice from VESTA; (b) Experimental XRD pattern for ZrB₂ [8]; (c) Lattice structure created and used for simulations.

3. Results and Discussion

3.1. Thermal Expansion Coefficient of ZrB₂

At near room temperatures (300 K), the CTE for all configurations coincide at 2.5~3.0 × 10⁻⁶ K⁻¹. Thermal expansion shows early stability between temperature ranges from 300 to 1200 K and suffers contraction at temperature exceeding 2000 K, approaching melting values. CTE for Hf impurity stands out in its rate of expansion and with notable breakdown at temperatures beyond 3000 K for all structures. NPT allows for dimension stabilization at given temperature which helps to predetermine the timesteps duration for thermodynamic equilibrium. Correlation showed early dependency on temperature which is consistent with reported analysis by Okamoto et al. [19]. However, with temperatures reaching ~1500 K, almost twice the Debye value for ZrB₂ (750 K), the expansivity peaked for all configurations (see Figure 4(a)). Hafnium (Hf) impurities in ZrB₂ was highest 7.9 × 10⁻⁶ K⁻¹ at ~1500 K while carbon (C), boron vacancies, boron interstitials defects were 6.1 × 10⁻⁵, 6.25 × 10⁻⁶, 7.7 × 10⁻⁶ K⁻¹ respectively. There is a slight correlation of CTE values to atomic radius (Hf > B > C) of the impurities added. The control simulation of perfect ZrB2 reported values within the ranges previously reported in experiments (6.66 - 6.93 × 10⁻⁶ K⁻¹) compared to simulations (5.9~6.68 × 10⁻⁶ K⁻¹) [8] [10] [12] [19]. Figure 4(a) also shows at temperatures with melting range (2900˚C - 3000˚C), there is a gradual breakdown, potentially structurally with uncontrollable expansivity. Also, the independent axial contribution to the expansion coefficient (see Figure 4(b)) was investigated using the normalized values of the lattice parameters. Figure 4(b) shows the normalized lattice parameter with very uniform slow contraction along x and overlapping dimensions for all crystal lattice configurations. Along z, one can attribute expansion and dimension offsets to structural changes. It aligns with the results obtained from X-ray diffraction study on ZrB₂ anisotropy with matching reflections [8].

3.2. Specific Heat Capacity and Structure of ZrB₂

To understand the effect of the impurities on the structural characteristics, the calculated total density of states (DOS) are plotted in Figure 5. The number of states for the perfect ZrB₂ reduced with increasing energy level. DOS for the covalently bonding impurities are closely overlapped with the pure crystal but with more dominating states likely from the projected orbitals of silicon and carbon at higher energy. The metallic impurities show states dominating at lower energy and slow increased. Potential bonding interaction between the covalent nature of ZrB₂ with both Si and C is noticeable in the overlap which is not the same for metallic impurities. Si extended the conduction band to 34.1 eV, while C impurities made negligible changes to the conduction band. Both Hf and W impurities extended both the conduction and the valence band. The conduction band in all cases are dominated by B orbitals. Increased width of valence band denotes increase in electron delocalization and therefore reduced band gap. On the other

Figure 4. (a) Coefficient of thermal expansion of ZrB₂ w/o defects; (b) Temperature dependency of lattice parameter w/o defects.

hand, decreasing width of conduction band indicates electron delocalization weakening. Occurrence of these impurities generated additional localized states and potential changes to conductive properties. Lattice contribution to heat capacity is computed from it’s the phonon DOS as shown below:

Figure 5. Phonon total density of states and partial density of states (a) for perfect ZrB₂ (b) for ZrB₂ with impurities of C, Si, Hf and W.

$C_v\left(t\right)=K_B{{\displaystyle \int }}^{\text{}}\frac{{\left(\hslash \omega /K_{B}T\right)}^2{\text{e}}^{\hslash \omega /K_{B}T}}{{\left[{\text{e}}^{\hslash \omega /K_{B}T}-1\right]}^2}F\left(\omega \right)\text{d}\omega$ (7)

where $F\left(\omega \right)$ represents the phonon DOS, lattice heat capacity ( $C_v$ ), $\hslash$ represents Planck’s constant and $k_B$ denotes Boltzmann’s constant.

The phonon dispersion for ZrB₂ structure was calculated along the high symmetry directions Γ-K-M-Γ as presented in Figure 6(a). The result for perfect crystal, shows dynamic stability as there are no imaginary modes, thus consistent with M. S. Daw et al. data [20]. Such stability is expected for structures with high atomic defect energy, and low tendency for atomic distortion. In the presence of impurities, there are observable differences. C and Si impurities show dynamical

Figure 6. Phonon dispersion (a) normal mode frequencies of perfect ZrB₂ (b) imaginary (dynamic) and normal modes for impurities.

stability similar to perfect structure but in addition to the dominant acoustic modes, there exist more low frequency optic modes. The added optic modes for Si are at relatively higher frequency than C. The case is different for the metallic impurities (Hf and W). As shown in Figure 6(b), negative phonons (imaginary branches) demonstrating dynamic instability for both configurations. It is noteworthy that some optic modes are also created by both Hf and W. The evolution of near zero frequency phonons transitioning, both normal and imaginary modes, points to changes in charge density area indicating that both of these metallic impurities have significantly distorted the original charge density distribution of ZrB₂. No phonon softening is exhibited in band structure due to increased weight from added impurities.

The computed specific heat is plotted for the respective impurities as show in Figure 7(a). The specific heat reaches a hard limit near the debye temperature (~750 K) following Dulong-Petit prediction. When lattice vibration reaches its maximum, phonon contribution will not increase even with temperature based on Debye model. Phonopy applies quasi-harmonic approximations (QHA) which is not applicable at very high temperatures of complete anharmonicity. Hence in this study, specific heat capacity in the temperature range of 0 - 1000 K was investigated. The presence of 5% vol of Hafnium (Hf) and Tungsten (W) impurities shows decreased specific heat capacity (158 J/mol/K) at reference 750 K temperature. This value is significantly lower than calculated for Carbon (C) and Silicon (Si) impurities. This can be described as stronger phonon scattering effects and reduced group velocity by Hf and W impurities in the ZrB₂ crystal. The subplot in Figure 7(a) also shows the offset between Si and C in C_v.

Figure 7(b) shows the changes in the thermodynamic properties with respect to lattice temperature. The dependences of free energy and entropy on lattice vibration and temperature are expressed in Equations 8 and 9. The internal energy changes in the material are clearly observed in entropy dispersed. It is non-trivial but shows differences in the microstates (degrees of freedom) of the various material. Translating it by equipartition theorem, the small offsets in heat capacity values (Figure 7(a)) is equivalent to changes in its degrees of freedom. The overlaps, further align with the differences in the phonon bands between the covalent (C and Si) and metallic impurities (Hf and W).

$F\left(T\right)=E_{tot}+E_{zp}+K_{B}T{\displaystyle \int F\left(\omega \right)\ln \left[1-{\text{e}}^{-\frac{\hslash \omega }{K_{B}T}}\right]\text{d}\omega }$ (8)

$S\left(T\right)=K_B\left\{{\displaystyle \int \frac{\frac{\hslash \omega }{K_{B}T}}{{\text{e}}^{\frac{\hslash \omega }{K_{B}T}}-1}F\left(\omega \right)\text{d}\omega }-{\displaystyle \int F\left(\omega \right)\left[1-{\text{e}}^{-\frac{\hslash \omega }{K_{B}T}}\right]\text{d}\omega }\right\}$ (9)

4. Conclusion

Having presented the approaches used for heat capacity and thermal expansion calculations with microstructural explanation using phonon properties, it has

Figure 7. (a) Overlay of calculated molar specific heat of ZrB₂ with impurities; (b) Thermodynamic properties: free energy, entropy and enthalpy of respective impurities with temperatures.

been shown that volumetric coefficient of thermal expansion is driven by anisotropic z-direction expansion with temperature. Impurity effects add to the knowledge gathering for better forecasting of thermal properties of ZrB₂-based materials. Exponential elevation of CTE near the melting temperatures shows structural breakdown and needs further investigation for phase transition. The challenge surrounding specific heat calculation, particularly using Phonopy is the quasi-harmonic level considerations, approximating anharmonic dependencies. Where available, computed thermal expansion and heat capacity data are compared to numerical and experimental data and are found to be within the range of experimentally reported results. Noteworthy is the redistribution of available states around the fermi level by metallic impurities, different from the covalent impurities in the electronic density of states. Further interest in structural differentiation in ZrB₂ and other transition metal diborides by analytical development is justified as the production process presents wide variations in density and purity level of the material.

Acknowledgements

This work used Oregon State University computing resources and also resources of the Bioinformatics Computational Research Group in University of California Riverside. Intel Corporation Education Assistance Program made this research work possible.

Data Availability Statement

The findings in this study are backed by computed data that are available from the corresponding author, [J.O. Ighere], upon reasonable request.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Rudy, E. (1988) Ternary Phase Equilibria in Transition Metal-Boron-Carbon Systems: Part V, Compendium of Phase Diagram Data. Technical Report AFML-TR-65-2. Wright Patterson Air Force Base (OH): Air Force Materials Laboratory.
[2]	Portnoi, K.I., Romasho, V.M. and Vyroshina, L.I. (1970) Phase Diagram of the Zirconium-Boron System. Poroshkoviaia Metallugia, 10, 68-71.
[3]	Mishra, S.K., Das, S. and Ramchandrarao, P. (2002) Microstructure Evolution during Sintering of Self-Propagating High-Temperature Synthesis Produced ZrB2 Powder. Journal of Materials Research, 17, 2809-2814. https://doi.org/10.1557/JMR.2002.0408
[4]	Gürcan, K. and Ayas, E. (2017) In-Situ Synthesis and Densification of HfB2 Ceramics by the Spark Plasma Sintering Technique Author Links Open Overlay Panel. Ceramics International, 43, 3547-3555. https://doi.org/10.1016/j.ceramint.2016.11.164
[5]	Zhang, S.C., Hilmas, G.E. and Fahrenholtz, W.G. (2006) Pressureless Densification of Zirconium Diboride with Boron Carbide Additions. Journal of the American Ceramic Society, 89, 1544-1550. https://doi.org/10.1111/j.1551-2916.2006.00949.x
[6]	Chubb, W. (1986) Coating a Uranium Dioxide Nuclear Fuel with a Zirconium Diboride Burnable Poison. US Patent Document 4,582,676/A, U.S. Commissioner of Patents, Washington DC.
[7]	Fiero, I.B. and Karoutas, Z.E. (2003) Implementation of Zirconium Diboride Burnable Absorber Coatings in CE Nuclear Power Fuel Assembly Designs. Westinghouse Non-Proprietary Class 3, WCAP-16072-NP.
[8]	Paxton, A.W., Ozdemir, T., Savkliyildiz, I., Bicer, H., Akdogan, K., Whalen, T., Zhong, Z. and Tsakalakos, T. (2016) Anisotrpic Thermal Expansion of Zirconium Diboride: An Energy-Dispersive X-Ray Diffraction Study. Journal of Ceramics, 2016, Article ID: 8346563.
[9]	Choi, H.J., Roundy, D., Sun, H., Cohen, M.L. and Louie, S.G. (2002) The Origin of the Anomalous Superconducting Properties of MgB2. Nature, 418, 758-760. https://doi.org/10.1038/nature00898
[10]	Zimmermann, J.W., Hilmas, G.E. and Fahrenholtz, W.G. (2008) Thermophysical Properties of ZrB2 and ZrB2-SiC Ceramics. Journal of the American Ceramic Society, 91, 1405-1411. https://doi.org/10.1111/j.1551-2916.2008.02268.x
[11]	Kaji, N., Skikano, H. and Tanaka, I. (1992) Development of ZrBr2-Graphite Protective Sleeve for Submerged Nozzle. Taikabustsu Overseas, 39-43.
[12]	Johnson, S.M. (2011) Ultra-High Temperature Ceramics: Application, Issues and Prospects. 2nd Ceramic Leadership Summit, Baltimore, 3 August 2011, 10.
[13]	Novikov, V.V., Matovnikov, A.V., Volkova, O.S. and Vasil’ev, A.N. (2017) Synthesis, Thermal and Magnetic Properties of RE-Diborides. Journal of Magnetism and Magnetic Materials, 428, 239-245. https://doi.org/10.1016/j.jmmm.2016.12.072
[14]	Mukubwa, A. (2018) Electronic Specific Heat of Iron Pnictides Based on Electron-Cooper Pair Interaction. Open Access Library Journal, 5, e5107.
[15]	Kresse, G. and Furthmuller, J. (1996) Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Physical Review B: Condensed Matter and Materials Physics, 54, 11169-11186. https://doi.org/10.1103/PhysRevB.54.11169
[16]	Kresse, G. and Joubert, D. (1999) From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Physical Review B: Condensed Matter and Materials Physics, 59, 1758. https://doi.org/10.1103/PhysRevB.59.1758
[17]	Plimpton, S.J. (1995) Fast Parallel Algorithms for Short-Range Molecular Dynamics. Journal of Computational Physics, 117, 1-19. https://doi.org/10.1006/jcph.1995.1039
[18]	Tersoff, J. (1988) Empirical Interatomic Potential for Silicon with Improved Elastic Properties. Physical Review B, 38, 9902. https://doi.org/10.1103/PhysRevB.38.9902
[19]	Okamoto, N.L., Kusakari, M., Tanaka, K., et al. (2003) Temperature Dependence of Thermal Expansion and Elastic Constants of Single Crystals of ZrB2 and the Suitability of ZrB2 as a Substrate for GaN Film. Journal of Applied Physics, 93, Article 88. https://doi.org/10.1063/1.1525404
[20]	Daw, M.S., Lawson, J.W. and Bauschlicher, C.W. (2011) Interatomic Potentials for Zirconium Diboride and Hafnium Diboride. Computational Material Science, 50, 2828-2835. https://doi.org/10.1016/j.commatsci.2011.04.038