Elastic stiffness tensors of Zr-$x$Nb alloy in presence of defects: A molecular dynamics study
Mohammad-Reza Basaadat Mahmoud Payami
cond-mat.mtrl-sci








School of Physics & Accelerators, Nuclear Science and Technology Research Institute, AEOI, 
P. O. Box 14395-836, Tehran, Iran

[mycorrespondingauthor]Corresponding author
mpayami@aeoi.org.ir


In a nuclear reactor, the Zr-xNb alloy, which is used as a structural material in the core region, is irradiated by energetic particles that cause the atoms to be displaced from their lattice sites and giving rise to crystal defects.  
The local changes in the atomic arrangements lead to local deformations of the solid and thereby changes of its local mechanical properties. Understanding the mechanisms behind this evolution in the core region of a reactor, and its monitoring or controlling is a critical task in nuclear industry. 
In this work, using extensive molecular dynamics 
simulations, we have studied the effects of radiation damage on the local mechanical properties of Zr-xNb alloy.    
In the first step, the effect of Nb-concentration on the mechanical stability of homogeneous 
Zr-xNb alloy is investigated. In the second step, we have studied the local changes of the elastic constants due to local changes of the microstructure. These local changes include presence and accumulation of vacancies in the form of dislocation loops or voids, accumulation of Nb atoms in the form of clusters of different morphologies. 
This study covers both cases of T=0^∘K and finite temperatures up to T=600^∘K.





Zr-xNb alloy; Crystal defect; Elastic constant; Vacancy; Dislocation loop; Self-interstitial; Molecular dynamics 








§ INTRODUCTION

Because of low cross-section for thermal-neutron capture, suitable corrosion-resistance in water, and good 
mechanical properties; zirconium and its alloys with niobium are widely used as structural materials in 
the core and fuel cladding of water-cooled nuclear 
reactors<cit.>. 

Pure zirconium is realized in three different crystal structures: the α-phase with an hcp lattice, 
being stable for temperatures lower than 860^∘ C; and the β-phase, having bcc crystal structure, 
being 
stable at higher temperatures<cit.>. It is also observed that the ω-phase, with 
hexagonal structure and space group p6/mmm, is obtained in a phase transition from α-phase at 
pressure of 2-7 GPa<cit.>. 
On the other hand, pure niobium crystallizes with a bcc structure at normal conditions. 

The Zr-xNb alloys at temperatures T<800^∘ K are characterized by α-phase crystalline 
structure<cit.>. 
Some of Zr-Nb based alloys that are widely used in nuclear reactors, are for example: Zr-1%Nb-O, known 
as M5, used in fuel cladding materials<cit.>; Zr-2.5%Nb-O alloys, used in CANDU 
reactors in Canada for pressure tube materials<cit.>; and a Zr-Nb-Sn-Fe alloy, known as 
ZIRLO, also used as a cladding material<cit.>.

Some prior experimental and theoretical investigations on the mechanical properties of Zr-Nb 
alloys have been performed by other researchers. For example: the elastic moduli of Zr single crystal and their temperature 
dependence were experimentally studied early 
in 1964 by Fisher and coworker<cit.>; the elastic constants of Nb single crystal were 
experimentally studied at 27^∘ C by Bolef<cit.>; the elastic constants of Nb-rich alloys 
Nb-xZr with bcc structure from 4.2^∘ K up to room temperature were experimentally determined by 
Hayes and coworker<cit.> in 1973; and also from zero temperature to the melting point were studied 
both experimentally and theoretically by Ashkenazi and coworkers<cit.> in 1978. 
Varshni<cit.>; 
Fast, et al.<cit.>; and Olsson<cit.> have investigated the temperature 
dependence of 
elastic constants of zirconium; Peng, et al. have studied 
the pressure effect on stabilities of self-interstitials in hcp-zirconium using first-principles  
method<cit.>; Liu, et al.<cit.>; Wang, et 
al.<cit.>; have calculated the elastic 
constants of niobium at normal pressure and high pressures, respectively; 
Kharchenko, et al. have investigated<cit.> the effect of Nb concentration on the 
lattice constants of Zr-xNb;  
Wang, et al.<cit.>; Weck, et al.<cit.>; Al-Zoubi, et 
al.<cit.>; have investigated the elastic constants 
for
Zr-Nb alloys; Xin, et al. have studied the point defect 
properties (such as formation energies) in hcp and bcc Zr with trace solute Nb by ab initio 
calculations<cit.>. However, to our knowledge, the mechanical properties of Zr-xNb 
alloys in presence of defects have not been thoroughly investigated, which will be studied in the present work.  

In the core of a nuclear reactor, Zr-Nb alloys are irradiated by energetic neutrons, high-energy photons, 
or energetic ions.
When the colliding particles are fast enough, they can transfer sufficient kinetic energies to the nuclei 
of atoms of the alloy. Due to this energy transfer, the atoms are ejected from their crystal lattice 
sites, leaving behind vacancies (V), diffuse throughout the material, and will finally stop within the 
crystal as self-interstitials (SI)<cit.>. The displaced atom with its corresponding 
vacancy forms a so-called Frenkel pair. 

At high enough temperatures and under the irradiation, the self-interstitial atoms and vacancies separately can diffuse throughout 
the crystal with different mobilities; 
giving rise to interstitial-vacancy recombinations or formation of larger atomic aggregates and vacancy 
clusters in the forms of dislocation loops or voids<cit.>. 
In addition, at intermediate temperatures the Nb atoms of the alloy may diffuse to a common region, 
leading to a local segregation of solute Nb atoms in a bcc structure<cit.>. 

The creation and evolution of the defects lead to local deformation of the alloy at those regions occupied by the 
defects. This deformation is due to local changes in the crystal lattice parameters or atomic arrangements 
which, in turn, lead to local changes in the elastic properties of the solid. 
Although the origins of such changes in the mechanical properties are known, the experimental 
determination of them in the core region of a reactor is a challenging task and therefore one must employ 
some theoretical models or resort to accurate simulation methods.

To determine the local changes in the mechanical properties of a defected crystal at different regions, we 
simulate that region in an MD box with a similar microstructure, and the size 
of the box is adopted sufficiently large so that the periodic images of the defects in simulation box do not 
interact.    

For large-scale molecular dynamics simulations, some interatomic potentials have already been developed 
and successfully used. For 
example: Mendelev, et. al.<cit.> using the embedded-atom method 
(EAM)<cit.> have developed 
Finnis-Sinclair type (“fs”) interatomic potential (Zr_3.eam.fs) for pure α-phase hcp 
zirconium, and 
have studied its phase transformations; Lin, et al.<cit.> using the EAM 
have developed an n-body potential for a Zr-Nb system and calculated some ground-state 
properties; and very recently, Smirnova, et al.<cit.> by adding angular-dependent 
term to the EAM 
potential, 
have developed a new 
“adp” interatomic potential (Zr_Nb.adp.txt) and have shown that the structure and properties of all 
Nb and Zr phases existing in the Zr-Nb alloy were reproduced with good accuracy. This potential is 
extensively used in the present study.

In this work, using the “adp” interatomic potential within the LAMMPS molecular-dynamics (MD) code package<cit.> we have performed extensive calculations on the structural and mechanical properties of  α-phase Zr-xNb systems with point and aggregated defects. 
The variations of the properties with respect to the temperature, up to working conditions of the reactor 
(600 ^∘ K) were also studied. 

The structure of this paper is as follows. Section <ref> is dedicated to the computational details; 
the calculation results are presented and discussed in section <ref>; and in section <ref>, we 
conclude this work. Finally, some useful relations employed for the calculations of mechanical properties are 
summarized in <ref>. 




§ COMPUTATIONAL DETAILS

To calculate the elastic constants of the solid with a given microstructure, 
we construct the simulation box with atomic configuration consistent with that microstructure and fully optimize the atomic positions as well as the cell geometry. 
In the next step, we take the relaxed supercell (simulation box) lattice vectors, denoted by {a_1 , a_2 , c}, and apply deformations according to each of the 6 deformation gradients {G_α} given by: 


    G_1= [ 1 + Δ    0     0 
     0     1     0 
     0     0     1 
    ],
    G_2= [ 1     0     0 
     0     1 + Δ    0 
     0     0     1 
    ],
    G_3= [ 1     0     0 
     0     1     0 
     0     0     1 + Δ
    ],


    G_4= [ 1    Δ    0 
    Δ    1     0 
     0     0     1 
    ],
    G_5= [ 1     0    Δ
     0     1     0 
    Δ    0     1 
    ],
    G_6= [ 1     0     0 
     0     1    Δ
     0    Δ    1 
    ],

and employ the Green-Lagrange strain tensor, ε^α defined by:

    ε^α=1/2(G^T_αG_α - I).

These 6 deformation gradients are applied one-by-one to the fully-relaxed structure such that only one independent deformation is considered at each time. For each of the 6 deformation modes, we perform positive and negative box displacements of "appropriate" magnitudes: Δ=-δ, +δ. For each of deformed structures, the stress tensor, σ^α, is calculated by MD run, allowing for relaxation of all the atomic degrees of freedom. Modifying the indices of strain and stress tensors according to Voigt-notation: 11 ↦ 1, 22 ↦ 2, 33 ↦ 3, 23 ↦ 4, 13 ↦ 5, 12 ↦ 6, the components of elastic constants C_ij^α were determined from the linear equation:


    σ_i^α =∑_j=1^6 C_ij^αϵ_j^α,     (i,α=1,2,⋯ ,6).


Finally, the nonzero components of elastic constants are determined from simple averaging over C_ij(-δ) and   
C_ij(+δ). In the hexagonal crystal system, it is common to take both lattice vectors a_1 and a_2 in the basal plane, whereas c orthogonal to the basal plane (See Fig. <ref>) so that we have C_11 = C_22. 





To run MD simulations at constant temperatures we utilize the Langevin thermostat<cit.> with sufficiently small damping factor.
Before applying deformations to the system, we first find the equilibrium box shape and size by employing the Parrinello-Rahman method in NPT (isothermal-isobaric) statistical ensemble imposing the constraints on pressure P=0, and off-diagonal components σ_ij=0; and then using the resulting optimized parameters we apply the deformations and continue our simulations employing NVT ensemble to determine the resulting stress components. Afterwards, using Eq.(<ref>), the components of elastic tensor are determined.  The equilibrations of such processes were achieved after 40 ps and 20 ps for NPT and NVT MD runs. The time-step was chosen as 0.001 ps.   

To determine the "appropriate" magnitude of box displacements in each case, we have tried the values δ=0.00001, 0.0001, 0.001, 0.01, 0.1; and plotted the resulting elastic constants (C_ij) as functions of δ,  and determined the appropriate range of values for which the results are insensitive to the choice of displacement magnitude, and then chose a value with low computation costs. For example, in Fig. <ref> we have shown the results for the case of Zr-1%Nb system for which δ=0.02 is a good choice.






Finally, the size of the simulation boxes were chosen to be sufficiently large to avoid the interactions between the specific defect in the box with those in the periodic image boxes. In addition, to simulate the properties at zero temperature, we have used a tiny value (near to zero) of T=0.01 K. So all listed properties at T=0 correspond to molecular dynamics simulations at this tiny value of temperature.



§ RESULTS AND DISCUSSIONS

To calculate the equilibrium geometry of pure α-Zr (i.e., Zr-xNb with x=0%) at T=0, a simulation box comprising of an 8× 8× 6 hcp unit cells containing 768 atoms was used. 
Then, allowing the changes of box size and shape as well as relaxing the atomic positions, we have fully optimized the geometry of the system and obtained the equilibrium lattice constants to be a=3.22Å and c=5.17Å with c/a=1.60, which are in good agreement with experimental values<cit.>: a=3.23Å and c=5.15Å with c/a=1.59.  In the second step, the elastic constants of Zr-xNb alloys with different Nb concentrations and morphologies are studied. In the next step, the effects of adding vacancies (V), self-interstitials (SI), vacancy dislocation loops, and voids on the elastic constants are investigated.    


 §.§ Homogeneous Nb distribution

As was mentioned earlier, in the course of time the Nb concentration of Zr-xNb may locally change and 
give rise to local changes 
of mechanical properties of the alloy. To see how the elastic constants depend on the local concentration x of Nb 
atoms, we have calculated the elastic constants of homogeneously distributed Nb atoms in 
Zr-xNb alloy for different values of x=0, 1, 2, 5, and 10%, where x=0 corresponds to pure Zr crystal. 

The variations of C_ij's with respect to Nb-concentration x at T=0 are shown in Fig. <ref>.



As is seen from Fig. <ref>, the rate of changes (in %) are 3, 17, 22, 57, and 325 for C_11, 
C_12, C_13, C_33, and C_44, respectively. Evidently, C_11 is the least sensitive 
elastic constant with 3% of change 
for 0≤ x≤ 10. The changes in other four elastic constants are relatively significant. C_12 shows 
decreasing behavior while other remaining three C_13, C_33, and C_44 increase with x. 

Using the conditions <ref>, it was verified that Zr-xNb alloys for x=1, 2, 5, and 10% are 
elastically stable. We therefore conclude that in a Zr-xNb alloy, there is no mechanical instability 
problems at inhomogeneous regions with 0≤ x≤ 10. However, it should be mentioned that although the 
relatively high concentration of Nb at a region in Zr-xNb alloy is mechanically harmless by itself, but 
this condition makes the region highly capable for hydrogen capture which decreases the elasticity and 
leads to development of cracks<cit.>. 

In Table <ref>, we have listed the Young's, bulk, and shear moduli (See <ref>) as given by Voigt 
and Reuss theories for polycrystalline material. Inspecting the listed elastic moduli, we observe that 
they all increase with increasing the x value in the alloy. This makes the alloy to become harder 
against any deformations. Increasing x by 10%, Young's modulus increases by 115% while the bulk and shear moduli increase by 12% and 126%, respectively.
 


Looking at the Poisson's ratio ν and anisotropy index A values listed in Table <ref>, we notice 
that both of them decrease by increasing the Nb concentration x in the alloy. That is, by 10% 
increase of x, the Poisson's ratio decreases by 15% from 0.429 to the value 0.364, while the “anisotropy 
index” decreases by 60% from 2.169 to the value 0.876.    

Because of thermal expansion of the alloy and thermal vibrations of atoms, the C_ij's undergo changes 
with 
increasing the temperature. They generally decrease with temperature, because the thermal expansion 
decreases the interacting forces between atoms. In Fig. <ref> we have shown the 
changes with temperature of elastic constants for 
homogeneous Zr-xNb including x=0, 1, 2, 5, and 10%. These changes are measured from the corresponding values at T=0.    



From Fig. <ref> we notice some general features that are worth mentioning: i)- The elastic 
constants C_11 and C_12 have similar behaviors in that for all studied concentrations x, they show decreasing behaviors with temperature, while for other elastic constants, the decreasing behavior is remained only for higher concentrations of x=5, 10%; ii)- C_13, 
C_33, and C_44 increase with temperature or remain almost constant for lower Nb-concentrations x=0, 1, 2%. In other words, the increasing behaviors of C_13, C_33, and C_44 for lower x values changes to decreasing behaviors for higher Nb concentrations.
A more careful inspection of Fig. <ref> reveals that except for C_12, other elastic constants 
have highest variations in the interval 0<T<600^∘ K for x=10%. In other words, the highest 
changes correspond to x=10% for C_11, C_13, C_33, C_44 by respective values of 11%, 
24%, 28%, 
32%.  



 §.§ Inhomogeneous Nb distribution

The local inhomogeneity of Zr-xNb alloy, produced by Nb aggregates, may be realized in different 
morphologies. To show how the morphologies can affect the local mechanical properties, we have considered 
three examples of different local inhomogeneities: planar Nb-aggregate (parallel to basal plane), three-dimensional hcp Nb-aggregate, and three-dimensional bcc Nb-aggregate (See Fig. <ref>).  



As shown in Fig. <ref>(a), to create a planar Nb-aggregate, we have replaced one layer (in 
basal plane) of 
Zr atoms in the simulation box with one layer of Nb atoms.
To create a 54-atom hcp-aggregate, as in Fig. <ref>(b) we have just replaced the Zr atoms in a 3× 3× 3 supercell with Nb atoms.  However, in Fig. <ref>(c) the replaced Nb-aggregate has bcc structure as is realized in precipitates.
The elastic constants obtained from the MD runs at T=0 for these three examples of Nb-inhomogeneities 
are presented in Table <ref>.  

In Table <ref>, the results show that C_11 decrease as we go from planar to hcp to bcc, whereas C_12 and C_13 decrease from planar to hcp but increase from hcp to bcc. Conversely, for C_33 and C_44 when we go from planar to hcp, they increase while by going from hcp to bcc they show decreasing behaviors. Quantitatively, the morphology of the Nb aggregates have contribution in the changes of elastic constants by at most ∼15GPa which corresponds to C_11 and C_33 whereas the least sensitive one is C_13 which varies at most by ∼2GPa.





In Fig. <ref>, we have plotted the temperature dependence of elastic constants for the systems shown in Fig. <ref>. The reference points C_ij(0) were taken as the corresponding values at T=0 (as tabulated in Table <ref>). The variations with respect to temperature are at most ∼20GPa which for C_11 belongs to planar and hcp inhomogeneities while for C_33 it belongs to local bcc Nb-precipitate.     



 §.§ Vacancies, self-interstitials, vacancy dislocation loops, and voids

Energetic particles when colliding the nuclei of atoms in a solid, are capable to eject atoms from their 
crystal lattice sites, leaving behind vacancies. On the other hand, the ejected atom move throughout the crystal and may stop somewhere as an interstitial atom [See Fig. <ref>(a)].    
The vacancies can diffuse throughout the crystal to join and give rise to formation of large vacancy 
clusters, also called voids [(See Fig. <ref>(b)] or in the form of vacancy dislocation loops as shown in Fig. <ref>(c). 





The existing internal stress fields of some extended defects may cause the single-atom vacancies and SI's or their clusters move in such a way that lead to partial segregation of atoms. However, in our study we do not consider their migration but have first investigated how the homogeneous distribution of V's or SI's can affect the local values of elastic constants.



 

To simulate the Zr-1%Nb alloy with 1%-V, we have removed homogeneously 8 Zr-atoms from the supercell containing 768 atoms, and for the 1%-SI, we have added 8 Zr-atoms homogeneously in the octahedral sites, so that the former and latter supercells have 768 and 776 atoms, respectively. Then the systems were geometrically optimized and the corresponding measured lattice constants are listed in Table <ref>.

From the results listed in Table <ref>, we see that adding 1%-V to the Zr-1%Nb system reduces (on the average) the values of a and c by 0.003 and 0.008Å, respectively; while adding 1%-SI increases them by 0.002 and 0.057Å, respectively. From these results we expect the interaction between atoms in the SI case to reduce compared to the other two cases, which in turn reduces the corresponding elastic constants. To verify this fact, we have obtained and compared the elastic constants for the three cases in Fig. <ref>.   
As is seen from Fig. <ref>, by adding the point defects, all C_ij's except for C_44 were decreased or not changed with respect to the perfect case, while C_44 increases in the case of V defect. In addition, as expected before, adding SI point defects to the system lead to higher changes in the C_ij values. 

Finally, as was mentioned, the vacancy point defects may migrate under the influence of stress fields of extended defects and join together to form vacancy clusters in the forms of dislocation loops or voids. To understand how these defects may affect the local mechanical properties, for a given concentration of vacancies, we have considered four arrangements of vacancy defects for a Zr-1%Nb system. In the first case, we have considered the single-atom vacancies (with 1% concentration) that are homogeneously distributed throughout the Zr-1%Nb alloy; in the second case, the vacancies that 
are joined to give rise to a maximum-sized dislocation-loop in basal plane; in the third case, dislocation-loop in prismatic plane; and in the fourth case, in the form of void. The corresponding elastic constants are listed in Table <ref>. 

From Table <ref>, two points are worth mentioning: First, introduction of 1%-V into Zr-1%Nb changes the values of C_11, C_12, C_13, C_33, and C_44 at most by 4, 6, 3, 6, and 16%, respectively; Second, as expected, vacancy cluster in the form of ⟨ c ⟩-loop decreases the C_33 with respect to perfect alloy, while in the case of ⟨ a ⟩-loop vacancy, the C_11 is decreased. 



 §.§ Eisentropic elastic constants

In this study, we have calculated the “isothermal” elastic constants at finite temperatures. However, in experimental data usually the “isentropic” values are reported and therefore it is necessary to calculate the corresponding conversion values. For hcp materials, an approximate relation between these two quantities are given by<cit.>:


    C_ij^S(T)=C_ij^T(T)+Δ_ij,

where C_ij^S(T) and C_ij^T(T) are “eisentropic” and “isothermal” elastic constants, respectively, and 

    Δ_ij=Tλ_iλ_jρ c_ϵ

with 

    λ_1=λ_2=α_a(C_11^T+C_12^T)+α_cC_13^T  ,  λ_3=2α_aC_13^T+α_cC_33^T.


In Eq.(<ref>), α_a and α_c are the linear thermal expansion coefficients in the 
directions of “a” and “c” axes of an hcp crystal, respectively. ρ and c_ϵ are the 
density and the specific heat at constant strain, respectively. In the α-phase Zr-xNb alloy we 
assume that the corresponding symmetries for an hcp perfect crystal are still valid and using the above 
relations we estimate the conversion values. The experimental value of c_ϵ for pure zirconium is 
reported<cit.> to be arround 270 Jkg^-1K^-1, and the linear thermal expansion 
coefficients α_a and α_c were experimentally reported<cit.> to be of order 
10^-6 K^-1. From our calculations, in the range of T=0^∘K and T=600^∘K, the average specific heat c̅_ϵ of pure α-Zr at P=0 is obtained to be 273 Jkg^-1K^-1; and for the linear expansion coefficients at the same temperature range, we have obtained the average values of α̅_a=2.3× 10^-6K^-1 and α̅_c=7.0× 10^-6K^-1. Putting these values in the 
conversion formulas, we estimate 
the value of order less than ∼ 1GPa which are practically insignificant. In fact, the conversion 
values would be significant for materials with high elastic constants, high thermal expansion 
coefficients, low densities, and low specific heats. 



§ CONCLUSIONS

In this study, employing an extensive MD simulations, the effects of different crystal defects on the elastic stiffness coefficients of Zr-1%Nb alloy were investigated.      
The defects considered in this study include: the local changes of Nb concentration, Nb clusters with different morphologies, single-atom vacancies and self-interstitials, vacancy clusters in the form of dislocation loops and voids. 
The local changes of Nb concentration up to 10% showed that C_11, C_12, C_13, C_33, and C_44
undergo changes of at most 3, 17, 22, 57, and 325%, respectively and that C_11 was the least sensitive 
elastic constant. In addition, it was shown that Zr-xNb alloys for x=1, 2, 5, 10% were 
elastically stable, and on the other hand, increasing the x value in the alloy makes it become harder 
against any deformations.
Moreover, studying the temperature effects on Zr-xNb showed that in the interval 0<T<600^∘ K, except for C_12, the highest 
changes in elastic constants correspond to x=10% for C_11, C_13, C_33, C_44 by respective values of 11%, 24%, 28%, 32%.  
In the context of Nb aggregates, it was shown that at T=0^∘K, the morphology has contributions by at most ∼15GPa in the elastic constants, which corresponds to C_11 and C_33; while, the variations with respect to temperature were at most by ∼20GPa for C_11 in planar form as well as hcp Nb-cluster, and for C_33 in the bcc Nb-precipitate. 
As to single-atom vacancies and self-interstitials, we have shown that adding 1%-V to the Zr-1%Nb system, the values of a and c reduces by 0.003 and 0.008Å, respectively; while addition of 1%-SI increases them by 0.002 and 0.057Å, respectively. 
To understand the effects of vacancy clusters in the forms of dislocation loops or voids on the local mechanical properties, we had considered four arrangements of vacancy defects for a Zr-1%Nb system with 1%-V. In the first case, the homogeneously-distributed single-atom vacancies were considered; in the second case, the maximum-sized dislocation-loop in basal plane; in the third case, dislocation-loop in prismatic plane; and in the fourth case, the maximum-sized void were considered. The results showed that firstly, introducing of 1%-V into Zr-1%Nb changes the values of C_11, C_12, C_13, C_33, and C_44 at most by 4, 6, 3, 6, and 16%, respectively; secondly, the vacancy cluster in the form of ⟨ c ⟩-loop decreases the C_33 with respect to perfect alloy, while in the case of ⟨ a ⟩-loop vacancy, the value of C_11 is decreased. 
Finally, it was shown that the difference between the calculated isothermal elastic constants and the experimental isentropic values were of order less than ∼ 1GPa for Zr-xNb alloy, which were practically insignificant and can be ignored. 



§ ACKNOWLEDGEMENT
 
This work is part of research program in School of Physics and Accelerators, NSTRI, AEOI.  




§ DEFINITION OF QUANTITIES USED IN MECHANICAL PROPERTIES




 §.§ Elastic constants

Deforming a crystal with a given strain ϵ_i, some forces appear tending to bring it to the equilibrium configuration (Voigt notation is used)<cit.>. Stress tensor is used for the description of these forces. For small values of strain, the stress tensor becomes proportional to the strain so that we will have:



    σ_i =∑_j=1^6 C_ijϵ_j,     (i=1,2,⋯ ,6)


The inverse of the Eq.(<ref>) is given by:


    ϵ_i =∑_j=1^6 S_ijσ_j,     (i=1,2,⋯ ,6)



Here, C_ij and S_ij define the elastic and the compliance tensors . At zero temperature, we can write:



    C_ij= ∂σ_i ∂ϵ_j,    S_ij= ∂ϵ_i ∂σ_j


The nonzero components of the elastic constant tensor can be derived using group theory. 
For hexagonal lattice, we have five nonzero independent components:

    C_11=C_22,   C_12,   C_13=C_23,   C_33,   C_44=C_55,   C_66=1  2(C_11-C_12).




 §.§ Elastic moduli

Young's modulus E describes the tensile elasticity, i.e., the tendency of an object to deform along an 
axis when opposing forces are applied along that axis; and it is defined as the ratio of tensile stress to 
tensile strain. 
The shear modulus (also called modulus of rigidity) G describes an object's tendency to shear (i.e., 
undergo shape deformation at constant volume) when applied upon opposing forces; and it is defined as the 
ratio of shear stress to shear strain[Care must be taken not to confuse the symbol G used here with those G_α used for deformation gradients.]. 
The bulk modulus, denoted by K, describes volumetric elasticity, i.e., the tendency of a solid to deform 
in all directions when it is uniformly loaded in all directions; and is defined as the ratio of volumetric 
stress to volumetric strain. Bulk modulus is the inverse of compressibility and is the extension of 
Young's modulus to three dimensional case.

The bulk and rigidity moduli for a polycrystalline material, that is macroscopically isotropic, are given 
from Voigt and Reuss theories by<cit.>:

    K_V = 1 9(C_11+C_22+C_33)+2 9(C_12+C_23+C_13)


    G_V = 1 15(C_11+C_22+C_33)-1 15(C_12+C_23+C_13)+1 5(C_44+C_55+C_66).


    1  K_R = (S_11+S_22+S_33)+2(S_12+S_23+S_13)


    1  G_R = 4  15(S_11+S_22+S_33)-4  15(S_12+S_23+S_13)+3  15(S_44+S_55+S_66).


For any stress and strain, we have:

    K_R≤ K≤ K_V;     G_R≤ G≤ G_V.


Poisson's ratio ν is defined as the ratio of transverse contraction strain to longitudinal extension strain in the direction of a stretching force. The Poisson's ratio ν and Young's modulus E for each scheme are given by:

    ν=1 2[1-3G/3K+G]      1/E=1/3G+1/9K.

From Eq.(<ref>) it is clear that:

    E_R≤ E≤ E_V.




 §.§ Anisotropy index

The Chung-Buessem empirical measure for anisotropy of shear modulus is given by<cit.>:

    A^G=G_V-G_R G_V+G_R,

and the universal index of anisotropy is defined as<cit.>:

    A^U=5(G_V  G_R)+(K_V K_R)-6.




 §.§ Elastic stability conditions

The necessary and sufficient elastic stability conditions for hexagonal crystal are given by<cit.>:

    Hexagonal:     C_11>  | C_12|,    (C_11+C_12)C_33> 2C_13^2,   
     C_44>0,       C_66>0.




§ DECLARATION OF COMPETING INTERESTS

The authors declare no affiliations with or involvement in any organization or entity with any financial or non-financial interest in the subject matter or materials discussed in this research paper.



§ DATA AVAILABILITY 

The raw or processed data required to reproduce these results can be shared upon an email to corresponding author (mpayami@aeoi.org.ir).
 


