3.9. Theory of Elasticity

download:pdf

3.9.1. Definition of Elastic Constants

To determine the elastic constants of a crystal, a deformation of the unit cell is created by changing the Bravais lattice vectors \(R = (a,b,c)\) of the undisturbed unit cell to \(R' = (a',b',c')\) using an engineering strain matrix \(e\)

(1)\[\begin{split}R' = R \begin{pmatrix} 1+e_{xx} \quad & \frac{1}{2}e_{xy} \quad & \frac{1}{2}e_{xz} \\ \frac{1}{2}e_{yx} \quad & 1+e_{yy} \quad & \frac{1}{2}e_{yz} \\ \frac{1}{2}e_{zx} \quad & \frac{1}{2}e_{zy} \quad & 1+e_{zz} \\ \end{pmatrix}\end{split}\]

The deformation leads to a change of the total energy of the crystal

(2)\[U = \dfrac{E_{tot}-E_{0}}{V_{0}} = \dfrac{1}{2} \sum_{i=1}^{6} \sum_{j=1}^{6} C_{ij} e_{i} e_{j}\]

where \(E_{0}\) is the total energy of the unstrained lattice, \(V_{0}\) is the volume of the undistorted cell and the \(C_{ij}\) are the elements of the elastic constant matrix with a notation that follows standard convention.

Both \(i\) and \(j\) run from 1…6 in the sequence \(\lbrace xx,yy,zz,yz,xz,xy \rbrace\). The tensor of elasticity has 36 elements, the elastic constants, but maximally 21 of these are independent. The symmetry of the unit cell can reduce the number of independent elastic constants for the various symmetry systems as shown in the following table:

Table 48 Number of independent elastic constants for unit cells of different symmetry
CRYSTAL SYSTEM ELASTIC CONSTANTS
Triclinic 21
Monoclinic 13
Orthorhombic 9
Tetragonal 6 or 7
Rhombohedral 6 or 7
Hexagonal 5
Cubic 3

For the sake of simplicity, in the so-called Voigt notation the six directions \(\lbrace xx, yy, zz, yz, zx, xy \rbrace\) are abbreviated by \(\lbrace 1, 2, 3, 4, 5, 6 \rbrace\). The following table describes the meaning each of the 21 elastic coefficients.

Table 49 The nomenclature of the 21 elastic coefficients.
ELASTIC COEFFICIENT STRAIN DIRECTION STRESS DIRECTION TYPE OF STRAIN VS STRESS
\(C_{11}\) \(xx\) \(xx\) Uniaxial stress in \(x\)-direction due to uniaxial strain in \(x\)-direction
\(C_{12}\) \(xx\) \(yy\) Uniaxial stress in \(y\)-direction due to uniaxial strain in \(x\)-direction
\(C_{13}\) \(xx\) \(zz\) Uniaxial stress in \(z\)-direction due to uniaxial strain in \(x\)-direction
\(C_{14}\) \(xx\) \(yz\) Shear stress in the \(y\)-\(z\)-plane due to uniaxial strain in \(x\)-direction
\(C_{15}\) \(xx\) \(zx\) Shear stress in the \(x\)-\(z\)-plane due to uniaxial strain in \(x\)-direction
\(C_{16}\) \(xx\) \(xy\) Shear stress in the \(x\)-\(y\)-plane due to uniaxial strain in \(x\)-direction
\(C_{22}\) \(yy\) \(yy\) Uniaxial stress in \(y\)-direction due to uniaxial strain in \(y\)-direction
\(C_{23}\) \(yy\) \(zz\) Uniaxial stress in \(z\)-direction due to uniaxial strain in \(y\)-direction
\(C_{24}\) \(yy\) \(yz\) Shear stress in the \(y\)-\(z\)-plane due to uniaxial strain in \(y\)-direction
\(C_{25}\) \(yy\) \(zx\) Shear stress in the \(x\)-\(z\)-plane due to uniaxial strain in \(y\)-direction
\(C_{26}\) \(yy\) \(xy\) Shear stress in the \(x\)-\(y\)-plane due to uniaxial strain in \(y\)-direction
\(C_{33}\) \(zz\) \(zz\) Uniaxial stress in \(z\)-direction due to uniaxial strain in \(z\)-direction
\(C_{34}\) \(zz\) \(yz\) Shear stress in the \(y\)-\(z\)-plane due to uniaxial strain in \(z\)-direction
\(C_{35}\) \(zz\) \(zx\) Shear stress in the \(x\)-\(z\)-plane due to uniaxial strain in \(z\)-direction
\(C_{36}\) \(zz\) \(xy\) Shear stress in the \(x\)-\(y\)-plane due to uniaxial strain in \(z\)-direction
\(C_{44}\) \(yz\) \(yz\) Shear stress in the \(y\)-\(z\)-plane due to shear strain in the \(y\)-\(z\)-plane
\(C_{45}\) \(yz\) \(zx\) Shear stress in the \(z\)-\(x\)-plane due to shear strain in the \(y\)-\(z\)-plane
\(C_{46}\) \(yz\) \(xy\) Shear stress in the \(x\)-\(y\)-plane due to shear strain in the \(y\)-\(z\)-plane
\(C_{55}\) \(zx\) \(zx\) Shear stress in the \(z\)-\(x\)-plane due to shear strain in the \(z\)-\(x\)-plane
\(C_{56}\) \(zx\) \(xy\) Shear stress in the \(x\)-\(y\)-plane due to shear strain in the \(z\)-\(x\)-plane
\(C_{66}\) \(xy\) \(xy\) Shear stress in the \(x\)-\(y\)-plane due to shear strain in the \(x\)-\(y\)-plane

3.9.2. Calculation of Elastic Constants

The simplest case by far is the cubic system where there are only three independent constants, \(C_{11}\), \(C_{12}\), and \(C_{44}\). We use this case to illustrate the manner in which the stiffness matrix elements may be determined from strain fields of the form of Eq. (1). If the applied engineering strain is \(e_{xx} = e\) with all other \(e_{ij}\) equal to zero, the energy change is \(U=C_{11}\dfrac{e^{2}}{2}\). This allows a unique determination of \(C_{11}\). If \(e_{yz} = e_{zy} = \frac{e}{2}\), with all other strain components zero, then \(U = C_{44} \frac{e^{2}}{2}\) and we have an independent determination of \(C_{44}\). The bulk modulus, \(B\), is the response to a uniform compression so applying the strain field \(e_{xx} = e_{yy} = e_{zz} = e\) allows the computation of \(B\) via the relation \(U=\frac{Be^{2}}{2}\). Similarly, the shear modulus can be calculated by using the strain field \(e_{zz} = e;e_{xx}=e_{yy}=-\frac{e}{2}\), whereupon \(U=3C'\frac{e^{2}}{2}\). Finally, the off-diagonal stiffness matrix element \(C_{12}\) can be calculated using one or other of the relations

(3)\[B = \dfrac{1}{3}(C_{11}+2C_{12})\]

and

(4)\[C' = \dfrac{1}{2}(C_{11}-C_{12})\]

Using both of these relations provides a useful independent check on the accuracy of the computation. A symmetry-general formulation of the calculation of elastic constants from total energy calculations is given by Le Page & Saxe [1].

3.9.2.1. Computational Aspects

The current implementation in MedeA is based on stress [2]. VASP calculates the stress tensor, i.e. the derivatives of the total energy of the cell with respect to changes of the 6 lattice parameters (cell lengths and angles), using analytic expressions. With a single value for the amount of distortion for each of the necessary engineering strains, it is possible to obtain the elastic constants. For simplicity, the amount of distortion is called “strain” in the VASP interface. A good value for the strain is 2%. If this value is too small the changes in the stress may be too small for accurate tracking. If e is too large, anharmonicity of the energy may play a role. Several strains can be used to get more points for the fitting procedure involved in the calculation of the elastic constants, thus leading to a higher numerical stability.

Hint

The elastic constants can be calculated for the equilibrium structure (at zero pressure) of a crystal, which is the usual case, or for distorted lattices. To ensure that the equilibrium structure does correspond to zero strain, the initial optimization of the theoretical cell parameters must be made with extremely high accuracy. In particular, for metallic systems a fine k-mesh (k-spacing of about 0.2/\({\mathring{\mathrm{A}}}\)) and the tetrahedron integration scheme should be used. The geometry convergence criterion should be 0.01 eV/\({\mathring{\mathrm{A}}}\) and the SCF convergence criterion needs to be set at 1.0E-06 or 1.0E-07.

Experimental values for the lattice constants, even where these are available, should not be used because the theoretical lattice will be in a state of strain for these values. This is because in the first instance, theory, and experiment will not agree exactly and in the second, the experimental error is unknown and probably of the same order as the strain field one needs to apply.

3.9.3. Polycrystalline Materials: Voigt-Reuss-Hill Method

The stiffness matrix calculated as above refers to the perfect lattice. However, approximate values referring to polycrystalline samples that take into account orientational disorder can be estimated from the single crystal data if one assumes that the elastic behavior is determined by the individual crystallites and not grain boundary effects. A thorough analysis of this issue can be found in a paper by Anderson [3]. The basic method, the Voigt-Reuss-Hill method, evolves from the expression for the shear modulus, \(G\),

(5)\[G_{Hill} = \dfrac{1}{2} \left[ G_{Voigt} + G_{Reuss} \right]\]

For cubic lattices, the Voigt and Reuss approximations for the shear modulus, \(G_{Voigt}\) and \(G_{Reuss}\), are:

(6)\[G_{Voigt} = \dfrac{1}{5}(C_{11} - C_{12} + 3C_{44})\]

and

(7)\[G_{Reuss} = \dfrac{5}{4(S_{11} - S_{12}) + 3S_{44}}\]

where \(S\) is the compliance matrix (i.e. the inverse matrix to \(C\)) and, therefore, the arithmetic mean should give a better estimate of \(G\) than either approximation. Checks by Anderson showed this to be true. Since the bulk moduli of the perfect crystal and a mono-phase polycrystal are the same, elastic data for the polycrystalline phase amount just to using the shear modulus as given by (5), (6), (7) in place of the value for the perfect crystal.

Once \(B\) and \(G\) are determined, Poisson’s ratio, \(\sigma\), and Young’s modulus, \(Y\), are given by the standard formulae of elasticity

(8)\[\sigma = \dfrac{1}{2} \dfrac{B-\frac{2}{3}G}{2B+\frac{1}{3}G}\]

and

(9)\[Y = \dfrac{9B}{1+\frac{3B}{G}}\]

With the Bulk and Young’s moduli, the Longitudinal modulus (or the P-wave modulus) can be estimated:

(10)\[M = \dfrac{3B(3B+Y)}{9B-Y}\]

Correspondingly, the mean values for the transverse and longitudinal sound velocities \(v_{t}\) and \(v_{l}\) are given in terms of the elastic moduli and the mean density \({\rho}\) by

\[v_{t} = \sqrt{\dfrac{G}{\rho}}\]

and

(11)\[v_{l} = \sqrt{\dfrac{B+\frac{4}{3}}{\rho}}\]

Finally, the mean sound velocity \(v_{m}\) is

(12)\[v_{m} = \left( \dfrac{1}{3} \left( \dfrac{2}{v_{t}^{3}} + \dfrac{1}{v_{l}^{3}} \right) \right)\]

3.9.4. Thermal Expansion and Heat Capacity

In the discussion in the previous sections, temperature played no role and the computation of the elastic constants referred to zero temperature. For a harmonic lattice the mean lattice parameters and elastic constants are independent of temperature. The observed temperature dependence of these quantities evolves from anharmonic behavior of the energy, which is why the computation of the coefficient of thermal expansion is one degree more difficult than the computation of the zero temperature elastic constants. The foundation for the computation of thermal corrections to lattice parameters was laid down by Grüneisen [4], remarkably, before the widespread acceptance of quantum mechanics! While at the time of Grüneisen even the simplest computation was beyond the capacity of hand calculators, today it is possible to implement Grüneisen’s theory, again using DFT as a generator of energies and forces, and the results can be very spectacular.

The full implementation of Grüneisen’s theory requires evaluation of the material’s phonon spectra and the derivatives of the individual phonon frequencies with respect to the cell parameters. However, there is a lower level of theory, the Grüneisen-Debye approximation, that allows computation of the thermal expansion coefficient with essentially no more effort than the elastic constants. In fact, for a cubic material, the only additional information to the elastic constants that is required is an averaged Grüneisen constant, \(\gamma_{G}\). As a first step, the Grüneisen constant could be set to a fixed value, e.g. to 2. The next step, could be the determination of \(\gamma_{G}\) from the pressure-volume, or internal energy-volume equation of state. A convenient approximation that incorporates \(\gamma_{G}\) as a parameter is [5]

(13)\[E(V) = \dfrac{BV_{0}}{\frac{5}{6} - \gamma_{G}} \left( \dfrac{V}{V_{0}} \right)^{\frac{5}{6} - \gamma_{G}} \left( ln \dfrac{V}{V_{0}} - \dfrac{1}{\frac{5}{6} - \gamma_{G}} \right) + E_{\inf}\]

where \(V_{0}\) is the equilibrium volume. Since \(\gamma_{G}\) is essentially the third derivative of \(E(V)\), the computation of this function must be carried out with high accuracy and great care, which is feasible.

Once \(\gamma_{G}\) is obtained, the only other input required to compute the thermal expansion coefficient is the Debye temperature, which is determined by exploiting its relationship to the sound velocity i.e.

(14)\[\theta_{D} = \dfrac{\hbar}{k_{B}} \left( \dfrac{6 \pi^{2} q}{V_{0}} \right) ^{\frac{1}{3}} v_{m}\]

where \(q\) is the number of atoms in the unit cell, \(V_{0}\) its volume, and \(\hbar\) and \(k_{B}\) are the Planck and Boltzmann constants, respectively. The lattice contribution to the specific heat capacity, \(C_{V}\), as a function of temperature, \(T\), can then be calculated via (see, for example, Mayer [5])

(15)\[C_{V}(T) = 9qk_{B} \left( \dfrac{T}{\theta_{D}} \right)^{3} \int_{0}^{x_{D}} x^{4} \dfrac{e^{x}}{(e{x} - 1)^{2}} dx\]

where \(X_{D} = \dfrac{\theta_{D}}{T}\) and the linear thermal expansion coefficient \(\alpha_{L}\) is given via the standard relation [6]

(16)\[\alpha_{L}(T) = \dfrac{1}{3} \gamma_{G} \dfrac{C_{V}(T)}{BV_{0}}\]

3.9.5. Isotropic Systems

In the isotropic regime, there are only two independent moduli and ,e.g., the Young modulus \(E\) can be expressed from Bulk \(K\), Shear \(G\) modulus, Poisson ratio \(\nu\), Lamé’s first coefficient \(\lambda\) or P-wave modulus \(M\) with the various relationships from below:

../../_images/image1421.png

The relationship between moduli and the elastic and compliance matrices is simple:

(17)\[C_{11} = M, \quad C_{12} = \lambda, \quad C_{44} = G \quad \textrm{and} \quad S_{11} = \dfrac{1}{E}, \quad S_{12} = \dfrac{\nu}{E}, \quad S_{44} = \dfrac{1}{G}\]

3.9.6. Elastic Constants for Disordered Cubic Systems

When reporting elastic data for a disordered structure, one can take inspiration from the averaging by Voigt and Reuss and “effective” coefficients for a simplified, anisotropic cubic system:

(18)\[\begin{split}\begin{matrix} \overline{C_{11}} = \dfrac{C_{11}+C_{22}+C_{33}}{3} \quad & \overline{S_{11}} = \dfrac{S_{11}+S_{22}+S_{33}}{3} \\ \overline{C_{12}} = \dfrac{C_{23}+C_{13}+C_{12}}{3} \quad & \overline{S_{12}} = \dfrac{S_{23}+S_{13}+S_{12}}{3} \\ \overline{C_{44}} = \dfrac{C_{44}+C_{55}+C_{66}}{3} \quad & \overline{S_{44}} = \dfrac{S_{44}+S_{55}+S_{66}}{3} \\ \end{matrix}\end{split}\]

This averaging is justified, as we can recreate the two matrices (in the simple anisotropic, cubic case) from the polycrystalline moduli and the Zener anisotropy ratio \(A\), which is independent of whether Voigt or Reuss averaging is chosen.

(19)\[\begin{split}\begin{matrix} \quad & Voigt \quad & & Reuss \\ K_{V} = & \dfrac{\overline{C_{11}}+\overline{2C_{12}}}{3} \quad & K_{R} = & \dfrac{1}{3\overline{s_{11}}+6\overline{s_{12}}} \\ G_{V} = & \dfrac{\overline{C_{11}}-\overline{C_{12}}+3\overline{C_{44}}}{5} \quad & G_{R} = & \dfrac{5}{4\overline{s_{11}} -4\overline{s_{12}} +3\overline{s_{44}}} \\ \end{matrix}\end{split}\]

And back from moduli to elastic matrix and compliance

(20)\[\begin{split}\begin{matrix} C_{11} = & K + \dfrac{10}{3G_{V}} - \dfrac{10G_{V}}{2A+3} & S_{11} = & \dfrac{1}{9K_{R}} + \dfrac{6}{6G_{R}} - \dfrac{5A}{2G_{R}(2+3A)} \\ C_{12} = & K - \dfrac{5}{3G_{V}} - \dfrac{5G_{V}}{2A+3} & S_{12} = & \dfrac{1}{9K_{R}} - \dfrac{5}{12G_{R}} + \dfrac{5A}{4G_{R}(2+3A)} \\ C_{44} = & \dfrac{5G_{V}}{2A+3} & S_{44} = & \dfrac{5A}{G_{R}(2+3A)} \\ \end{matrix}\end{split}\]

3.9.7. Hashin-Shtrikman bounds for cubic systems

Voigt [7] (1928) averaged the elastic stiffnesses (\(C_{ij}\)) over all space. Reuss [8] (1929) averaged the elastic compliances (\(S_{ij}\)). Hill [9] (1952) showed that these averaged values are the least upper bound and the greatest lower bound, respectively, for the aggregate.

The Voigt and Reuss averages are given by:

(21)\[ \begin{align}\begin{aligned}\begin{matrix}\\\begin{split} & Voigt & \quad \quad & Reuss \\ \overline{C}_{11} = & \dfrac{C_{11}+C_{22}+C_{33}}{3} & \quad \quad \overline{S}_{11} = & \dfrac{S_{11}+S_{22}+S_{33}}{3} \\\end{split}\\\begin{split} \overline{C}_{12} = & \dfrac{C_{23}+C_{13}+C_{12}}{3} & \quad \quad \overline{S}_{12} = & \dfrac{S_{23}+S_{13}+S_{12}}{3} \\\end{split}\\\begin{split} \overline{C}_{44} = & \dfrac{C_{44}+C_{55}+C_{66}}{3} & \quad \quad \overline{S}_{44} = & \dfrac{S_{44}+S_{55}+S_{66}}{3} \\\end{split}\\\begin{split} K_{V} = & \dfrac{\overline{C}_{11}+2\overline{C}_{12}}{3} & \quad \quad K_{R} = & \dfrac{1}{3\overline{S}_{11}+6\overline{S}_{12}} \\\end{split}\\\begin{split} G_{V} = & \dfrac{\overline{C}_{11}-\overline{C}_{12}+3\overline{C}_{44}}{5} & \quad \quad G_{R} = & \dfrac{5}{4\overline{S}_{11}-4\overline{S}_{12}+3\overline{S}_{44}} \\\end{split}\\ \end{matrix}\end{aligned}\end{align} \]

Hashin and Shtrikman [10] showed that the Voigt and Reuss bounds could be improved, and they developed expressions for the corresponding bounds of aggregates of cubic crystals. As it turns out, the difference between the upper and lower bounds for cubic crystals is improved by at least an order of magnitude for most materials. For a single-phase aggregate of a cubic material, the bulk modulus, \(K\), is given unambiguously by

(22)\[K = \dfrac{1}{3}(\overline{C_{1}}+2\overline{C_{12}})\]

And the rigidity is bounded by

(23)\[ \begin{align}\begin{aligned}\begin{matrix}\\\begin{split} G_{1}^{\ast} & = & G_{1} + \dfrac{3}{\frac{5}{G_{2}-G_{1}}-4\beta_{1}} \quad & \quad G_{2}^{\ast} & = & G_{2} + \dfrac{2}{\frac{5}{G_{1}-G_{2}}-6\beta_{2}} \quad \\\end{split}\\\begin{split} G_{1} & = & \dfrac{1}{2}(C_{11}-C_{12}) \quad & \quad \quad \quad G_{2} & = & C_{44} \quad \\\end{split}\\\begin{split} \beta_{1} & = & - \dfrac{3(K+2G_{1})}{5G_{1}(3K+4G_{1})} \quad & \quad \beta_{2} & = & - \dfrac{3(K+2G_{2})}{5G_{2}(3K+4G_{2})} \quad \\\end{split}\\\begin{split} \quad & = & - \dfrac{2}{15} \dfrac{4C_{11}-C_{12}}{C_{11}(C_{11}-C_{12})} \quad & \quad & = & - \dfrac{C_{11}+2C_{12}+6C_{44}}{5C_{44}(C_{11}+2C_{12}+4C_{44})} \quad \\\end{split}\\\end{matrix}\end{aligned}\end{align} \]

The factors \(G_{1}\) and \(G_{2}\) are no coincidence, remember that in isotropic elasticity, the shear modulus \(G = C_{44}\), which turns up here as \(G_{2}\). The same ingredients appear in Zener’s anisotropy coefficient \(A\), which can be written as

(24)\[A = \dfrac{G_{2}}{G_{1}}= \dfrac{2C_{44}}{C_{11}-C_{12}}\]

This also determines which of the two values is the upper and lower bound:

(25)\[\begin{split}\begin{matrix} & A<1: & G_{H}=G_{1}^{\ast} & G_{S}=G_{2}^{\ast} \\ & A=0: & G_{H}=G_{S}=G_{Hill} & \\ & A>1: & G_{H}=G_{2}^{\ast} & G_{S}=G_{1}^{\ast} \\ \end{matrix}\end{split}\]

The larger term of \(G_{1}^{\ast}\) and \(G_{2}^{\ast}\) is termed the Hashin rigidity, \(G_{S}\), the smaller Shtrikman rigidity, \(G_{H}\). These two are much closer to each other then Voigt and Reuss bounds in equation

(26)\[C_{V} = \dfrac{\overline{C}_{11} - 2 \overline{C}_{12} + 3 \overline{C}_{44}}{5} > G_{H} > G_{Hill} > G_{S} > G_{R} = \dfrac{5}{4 \overline{S}_{11} - 4 \overline{S}_{12} + 3 \overline{S}_{44}}\]

3.9.8. Hill-Walpole Bounds for Amorphous Systems

Atomistic simulations of amorphous glassy materials are typically performed using periodic cells with lateral dimensions of no more than a few tens of Angstroms. Moreover, to model accurately and to predict the expected macroscopic properties of a given material it will often be necessary to sample many distinct realizations (i.e. cells), each differing in local packing of atoms. In the case of mechanical properties, computed 6x6 matrices of stiffness or compliance constants may differ quite considerably between cells and each may be slightly anisotropic, even when the parent material is homogeneous and perfectly isotropic.

An important question confronting modelers seeking to predict mechanical properties and often to quantify differences between materials, concerns how to combine elastic constants obtained for nanoscopic regions represented by amorphous cells to define the macroscopic elastic constants as precisely as possible. Instinctively, the basic approach involves viewing the macroscopic material as a composite made up of many distinct domains - i.e. individual cells - embedded in a matrix formed by all other cells.

As well-established in composite theory, the extreme upper and lower bounds of the elastic constants can then be determined by computing averages of the individual stiffness and compliance matrices. These are the so-called Voigt and Reuss bounds, which are obtained by assuming that each individual domain is subjected to the same strain or the same stress respectively.

(27)\[\langle C \rangle _{Voigt} = \dfrac{1}{n} \sum_{i} C_{i} \quad \textrm{and} \quad \langle S \rangle _{Reuss} = \dfrac{1}{n} \sum_{i} S_{i}\]

Since in reality neither assumption is rigorously true, these represent extreme bounds, as the forces between domains in the Voigt average cannot be in equilibrium; likewise the distorted domains in the Reuss average cannot fit together. It is just the most simplistic estimate of the likely macroscopic properties, though it has been suggested that simply taking the mid-point of the bounds interval can sometimes provide useful estimates of the material elastic constants.

In the case of amorphous glassy materials as frequently modeled by a collection of cells containing a few thousand to a few tens of thousands of atoms, where lateral dimensions range from around 20 to 50 Angstroms, upper and lower bounds calculated in this way typically show large differences. With the goal of improving the reliability of predictions of the elastic constants of macroscopic samples of material, Suter and Eichinger [11] accordingly investigated the possibility of obtaining refined bounds estimates from raw ‘amorphous cell’ elastic constant data.

By again proceeding with the analogy with composite models, and taking the view that each domain can be considered as being embedded within a homogeneous continuum with elastic properties equal to the average material, these authors applied theory developed earlier by Hill [12] and by Walpole [13] to compute improved and typically much narrower bounds estimates.

Suter and Eichinger then proceeded to investigate the effectiveness of the approach at ‘predicting’ the elastic constants of amorphous glassy polystyrene using input data ‘synthesized’ by adding random fluctuations to actual experimental stiffness coefficients, where the random fluctuations were chosen based on those observed in typical atomistic simulations. The results essentially demonstrated that the range of the computed new ‘Hill-Walpole’ bounds was significantly smaller than the difference of the Voigt-Reuss averages, and consistently closer to the known values used to generate the input data than the midpoint of the Voigt-Reuss range.

[1]Y. Le Page and P. Saxe, Symmetry-General Least-Squares Extraction of Elastic Coefficients From Ab Initio Total Energy Calculations, Phys. Rev. B 63, 174103 (2001).
[2]Y. Le Page and P. Saxe, Symmetry-General Least-Squares Extraction of Elastic Data for Strained Materials From Ab Initio Calculations of Stress, Phys. Rev. B 65, 104104 (2002).
[3]O. L. Anderson, A Simplified Method for Calculating the Debye Temperature From Elastic Constants, J. Phys. Chem. Solids 24, 909 (1963).
[4]E. Grüneisen, Zustand des festen K&ouml;rpers, in: Thermische Eigenschaften der Stoffe, Handbuch Der Physik, (Springer Berlin, Heidelberg, 1928).
[5](1, 2) B Mayer, H Anton, E Bott, M Methfessel, J Sticht, J. Harris, and P. C. Schmidt, Ab-Initio Calculation of the Elastic Constants and Thermal Expansion Coefficients of Laves Phases, Intermetallics 11, 23 (2003).
[6]N. W. Ashcroft and N. D. Mermin, Solid State Physics, Holt-Saunders, 1976.
[7]W. Voigt, Lehrbuch Der Kristallphysik (Mit Ausschluss Der Kristalloptik), (Teubner, Leipzig, 1928).
[8]A. Reuss, Berechnung der Fliessgrenze von Mischkristallen auf Grund der Plastizit&auml;tsbedingung f&uuml;r Einkristalle, Z. Ang. Math. Mech. 9, 49 (1929).
[9]R. Hill, The Elastic Behaviour of a Crystalline Aggregate, Proc. Phys. Soc.. Sect. A 65, 349 (1952).
[10]Z. Hashin and S. Shtrikman, A Variational Approach to the Theory of the Elastic Behaviour of Polycrystals,” J. Mech. Phys. Solids 10, 343 (1962); On Some Variational Principles in Anisotropic and Nonhomogeneous Elasticity, J. Mech. Phys. Solids 10, 355 (1962).
[11]U. W. Suter and B. E. Eichinger, Estimating Elastic Constants by Averaging Over Simulated Structures, Polymer 43, 575 (2002).
[12]R: Hill, Continuum Micro-Mechanics of Elastoplastic Polycrystals, J. Mech. Phys. Solids 13, 89 (1965); A Self-Consistent Mechanics of Composite Materials, J. Mech. Phys. Solids 13, 213 (1965).
[13]L. J. Walpole, On Bounds for the Overall Elastic Moduli of Inhomogeneous Systems-I, J. Mech. Phys. Solids 14, 151 (1966); On the Overall Elastic Moduli of Composite Materials, J. Mech. Phys. Solids 17, 235 (1969).
download:pdf