A comparison of polycrystalline elastic constants computed by analytic homogenization schemes and FEM

BCC magnesium-lithium alloys are a promising light-weight structural material.  As a first step in a theoretically guided materials design strategy single crystal elastic constants for bcc magnesium-lithium alloys with different compositions were computed using ab initio methods.  These single crystal elastic constants were then used to predict the corresponding polycrystalline elastic constants using various analytic homogenization techniques Voigt, Reuss, and a self-consistent approach) as well as}{the finite element method.  As expected, the Voigt and Reuss bounds form the upper and lower bounds on the polycrystalline elastic constants, which the predicted values of the self-consistent approach and finite element fall in between.  Additionally, the difference between the elastic constants derived from the self-consistent approach and the finite element method is small illustrating the power and value of the self-consistent approach.

Introduction
Magnesium-lithium alloys have the potential to be one of the lightest possible metallic alloy systems.  The density of MgLi alloys is expected to range between 1.74 g/cm3 (Mg) and 0.58 g/ cm3 (Li).  An additional benefit of alloying Mg with Li is that Li stabilizes the body-centered-cubic (BCC) structure over the hexagonal closed packed (HCP) structure with as little as 30 at. % Li.  Since BCC materials are generally more ductile at room temperature than HCP materials, they are favored in manufacturing operations that require room or low-temperature forming operations.  
While there is little doubt that bcc MgLi alloys will be light-weight, it is unclear whether these alloys will offer advantageous polycrystalline elastic properties.  As part of a theoretically guided materials design strategy, the single crystal elastic constants of various BCC-MgLi alloys were calculated using density-functional theory (DFT) \cite{Art08}.  However, there remains the problem of calculating the macro-scale polycrystalline elastic properties (like Young's and shear modulus) from the single crystal elastic components (Cij's) calculated at the atom level.
It is possible to estimate the polycrystalline elastic constants using various analytic homogenization techniques (like Voigt, Reuss, and the self-consistent approach) as well as the finite element method (FEM).  These analytic approaches make assumptions on the stress-strain state or geometry of the polycrystal in order to derive simple closed form solutions.  The advantage of FEM is that no such assumptions are neseccary.  In this paper, the usefullness of the analyitc estimates of polycrystalline elastic constants for materials with cubic symmetry is investigated.  The polycrystalline Young's modulus, shear modulus, and Poisson's ratio of BCC MgLi alloys are estimated using both analytic homogenization techniques and FEM from the same single crystal data and then compared.

Computational Methods

Ab Initio Calculation of Single Crystal Elastic Constants
DFT calculations \cite{Hohenberg1964,Kohn1965} employing the generalized gradient approximation \cite{Perdew1996} were performed using a plane wave pseudopotential approach as implemented in the Vienna Ab-initio Simulation Package (VASP) code \cite{Kresse1993,Kresse1996}.  Binary alloys were described by cubic 2x2x2 supercells containing a total of 16 atoms. The plane wave cutoff energy was set to 260 eV and a 16x16x16 Monkhorst-Pack mesh was used to sample the Brillouin zone.  The ground state properties (0 K) of pure bcc Li and Mg as well as 11 alloy compositions were studied.  These alloys were created by systematically replacing Mg and Li atoms in such a way that preserved cubic symmetry. Alloy compositions ranged from 6.25 at. %  Li (one Li atom in a 16 atom supercell) to 93.75 at. % Li (15 Li atoms in a 16 atom supercell). 
The three independent elastic constants of a single-crystal in the (100) direction with cubic symmetry were determined by distorting the bcc supercell in three different ways \cite{Chen2003}.  The first distortion is the isotropic volume change from which the single crystal bulk modulus is determined.    The single crystal bulk modulus for materials with cubic symmetry can be expressed as a linear combination of C11 and C12

The second distortion into a tetragonal structure (tet) and the third distortion into a trigonal one (tri) are defined as

where $\delta$ is the distortion.  To avoid non-linear contributions, $\delta$ was limited to $\pm$4\%.  Both tet and  tri lead to changes in the total energy density of the system as a function of $\delta$.  This change in energy density can then be used to calculate the single crystal elastic constants via

where Utet is the strain energy density due to tet  and Utri is the strain energy density due to tri.  Five different tetragonal and five trigonal distortions were used to calculate the corresponding single crystal elastic constants.

Analytic Homogenization - Voigt, Reuss, and Self-consistent
The procedure to estimate polycrystalline elastic constants begins by first estimating the polycrystalline bulk modulus (B*) and shear modulus (G*) from single crystal values.  Assuming that all possible grain orientations are possible and equally likely (i.e. a non-textured material), both B* and G* can be used to calculate Young's modulus (Y*) and Poisson's ratio (*) for an elastically isotropic polycrystal via well-established elasticity relations.
The upper bound on B* and G* was first calculated by Voigt \cite{Voigt1928} based on the assumption that during elastic deformation a polycrystal adopts a local uniform strain state.  The lower bound on B* and G* was calculated by Reuss \cite{Reuss1929} based on the assumption that during elastic deformation the polycrystal adopts a local uniform stress state.  The uniform strain and stress assumptions used by Voigt and Reuss represent two extremes that are rarely observed. For materials with cubic symmetry, the Voigt and Reuss bounds on B* (B*V and B*R respectively) are identical meaning that the polycrystalline bulk modulus is equal to the single crystalline bulk modulus
The Voigt and Reuss bounds on the polycrystalline shear modulus (G*V and G*R respectively) for materials with cubic symmetry depend on the single crystal elastic components via

where Cij are terms in the stiffness matrix and Sij = Cij-1.
A more realistic estimate of G* can be made by employing a self-consistent scheme.  Hershey \cite{Hershey1954} derived a closed form solution to the self-consistent homogenization of the shear modulus (G*SC) for materials with cubic symmetry.  Considering a polycrystal made up of spherically shaped grains of equal size, Hershey showed that G*SC could be calculated by solving the following quartic equation

Once homogenized values of G* and B* have been estimated, other elastic constants of an elastically isotropic polycrystalline aggregate can also be calculated with standard elasticity relations.  In particular,Y* and * are

Finite Element Method Homogenization
Polycrystalline elastic constants were calculated using the anisotropic elastic material routine within the commercial finite element code MSC.Marc200x.  In order to investigate the effect of grain size distribution, two different polycrystals were meshed and simulated.  The first polycrystal contained 96 (4 grain x 6 grain x 4 grain) equally shaped square grains.  Each grain was meshed with 27 quadratic brick elements resulting in a mesh with 2592 total elements.  The second polycrystal contained 84 grains of varying sizes and shapes.  All the grains in this polycrystal were meshed with quadratic brick elements and the resulting mesh had a total of 4096 elements.  Both the square grain and non-uniform grain size polycrystal meshes are shown in Fig 1.  

In both polycrystals, each grain was assigned a random orientation and the elasticity tensor (Cijkl) containing the single crystal elastic components was then rotated into the corresponding grain's orientation (i.e. the elastic constants in each grain differed only by a rotation).

The FEM determined polycrystalline Young's modulus (Y*FEM) and Poisson's ratio (*FEM) were calculated by simulating a uniaxial tensile test.  A fixed displacement normal to the y-direction in the polycrystal was prescribed while displacements on three other orthogonal planes were fixed to prevent the mesh from translating.  The overall stress-strain response was then calculated from the nodal reaction forces and displacements.  Furthermore, because polycrystals containing around 100 grains constitutes a small stastitical population, stress-strain response of 5 different polycrystals was simulated and an average Y*FEM and *FEM was calculated.  The G*FEM was then calculated assuming elastic isotropy.

Results
Single crystal elastic constants for ordered bcc MgLi alloys with compositions ranging from 0-100 at. % Li were calculated using ab initio methods.  It should be noted that bcc MgLi alloys are thermodynamically stable for Li concentrations greater than 30 at. \%.  Nevertheless, the bcc elastic constants were predicted in regions where the bcc crystal structure is not stable in order to get the complete picture on their concentration dependence. The variation of the three independent elastic components (C11, C12, and C44) is shown in Fig 2}. 

While there is little literature data available to validate these single crystal elastic constants, it is possible to compare ab initio calculated and experimental single crystal elastic constants for pure bcc Li at 78 K \cite{Nash1959}.  The difference between the experimental and predicted elastic components is between 1-6 % (or approximately 1 GPa).
These single crystal elastic components were then used as the basis for all the analytic and FEM estimates of polycrystalline elastic constants.  The dependence of Y*, G*, and * as a function of Li content are shown in Figs.3-5.  
  

Fig. 5 -The dependence of polycrystalline Poisson's ratio on lithium content for bcc MgLi.

In each case, the Voigt and Reuss bounds form upper and lower bounds, as should be the case.  The difference between these two bounds varies between 20-130% for G* and Y* and 25-50% for *.  Such large differences between the upper and lower bounds illistrate that these bounds by themselves do not necessarily restrict the elastic constants to a small-sub set from which reasonble values can be extracted.  Intrestingly enough, the self-consistent and FEM predictions fall close to the midpoint between upper and lower bounds.  This result validates Hill's \cite{Hill1952} suggestion that B* and G* can be reasonably estimated by taking an arithmetic mean of the upper and lower bound.
In Figs. \ref{poly-Y} and \ref{poly-G}, the Voigt and Reuss approaches estimate polycrystalline elastic constants for bcc Mg, but the self-consistent of FEM approaches do not.  In the self-consistent case, all of the roots to Eq.~\ref{hershey-orig} are imaginary; while for FEM, the simulations crash immediately due to a non-positive definite stiffness matrix.  The inability of the self-consistent and FEM methods to calculate real elastic constants for bcc Mg indicates that it is mechanically unstable.  In other words, bcc Mg will undergo a phase transformation upon the application of strain.  The negative G*V and Y*V also hints at this mechanical instability.  These results are in agreement with the ab initio Bain path calculations by Jona and Marcus \cite{JonaMarcus2002}.
The elastic constants predicted by the self-consistent approach and do not differ greatly from those predicted by FEM.  FEM in general predicts a slightly lower modulus value and a slightly higher * value.  The difference between the predicted G* and Y* values from both approaches was between 2-17% and between 0-4% for *.  However, percent differences are not be the best way to compare the two approaches because the difference between the self-consistent and FEM values is constantly between 1-3 GPa for shear and Young's modulus and between 0.003-0.02 for Poisson's ratio.  Therefore, alloys with low G* and Y* values (like those with high Li content) will have a higher % difference than those alloys with higher values of G* and Y*.  This constant error trend indicates that FEM should be used carefully when investigating the elastic properties of soft polycrystalline materials.  Despite this drawback, the ability to predict similar polycrystalline elastic constants (within 1-3 GPa) using a quartic equation (solved instantly) and FEM (solved on the order of minutes to hours) illustrates the power and value of the self-consistent approach.
A comparison of the polycrystalline elastic constants predicted by FEM using a uniform and non-uniform grain size distribution with those predicted using the self-consistent approach reveals that using a non-uniform grain size distribution results in elastic constants that are closer to the self-consistent predictions.  The differences between G*, Y*, and * for both FEM approaches and the self-consistent approach are summarized in Table \ref{tab:fem}.  In general, the FEM predictions using a non-uniform grain size polycrystal are 50% closer to the self-constistent results than the FEM results based on a polycrystal containing a uniform grain size.  In one sence, this overall trend was unexpected since the self-consistent approach was dervied from a polycrystal with a uniform spherical grain size distribution.  However, FEM simulations based on a more realistic polycrystal should lead to a closer prediction of the exact solution indicating that the exact solution lies very close to values predicted using the self-consistent approach.

Conclusions

Single crystal elastic components (C11, C12, and C44) were predicted for 11 bcc MgLi alloys as well as bcc Mg and Li using ab initio methods.  These elastic components were then used to estimate polycrystalline elastic constants (Y*, G*, and *) using both analytic homogenization and FEM methods.  As expected, the upper and lower bounds for each of the polycrystalline elastic constants was defined by the Voigt and Reuss solutions.  While the difference between these two bounds was too large to be useful, the arithmetic mean of the two bounds did provide a reasonable approximation.  
For elastically isotropic materials with cubic symmetry, the difference between the polycrystalline elastic constants predicted using and FEM and those predicted using an analytic self-consistent approach is small.  Y* and G* computed using FEM were 1-3 GPa smaller than those computed using the analytic self-consistent approach, while the FEM based values for * were 0.00-0.15 higher.  The effect of grain size distribution on the FEM results was also investigated.  The FEM predictions using a non-uniform grain size polycrystal are 50% closer to the self-consistent results than the FEM results based on a polycrystal containing a uniform grain size.  The good agreement between the polycrystalline elastic constants derived from the self-consistent approach and FEM illustrates the power and value of applying the self-consistent approach to materials with cubic symmetry.

References
[1] W. A. Counts, M. Fri´ak, D. Raabe, J. Neugebauer, Submitted
to Acta Materilia.
[2] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
[3] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[4] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett.
77, 3865 (1996).
[5] G. Kresse and J. Hafner, Phys. Rev. B 47(1), 558–561
(1993).
[6] G. Kresse and J. Furthm¨uller, Phys. Rev. B 54(16), 11169–
11186 (1996).
[7] K. Chen, L.R. Zhao, and J. S. Tse, J. Appl. Phys. 93(3),
2414–2417 (2003).
[8] W. Voigt, Lehrbuch der Kristallphysik (Teubner, Stuttgart,
1928).
[9] A. Reuss, Z. Angew. Math. Mech. 9, 49 (1929).
[10] A.V. Hershey, J. Appl. mech. 9, 49 (1954).
[11] T. Massalski, Binary Alloy Phase Diagrams (ASM International,
1990).
[12] H. Nash and C. Smith, J. Phys. Chem. Solids 9, 113 (1959).
[13] R. Hill, Proc. Phys. Soc. Lond. A 65(349) (1952).
[14] F. Jona and P.M.Marcus, Phys. Rev. B 66(094104) (2002).