Elastic Properties

The Bulk Modulus

The bulk elastic properties of a material determine how much it will compress under a given amount of external pressure. The ratio of the change in pressure to the fractional volume compression is called the bulk modulus of the material.

A representative value for the bulk modulus for steel is 160 x 10⁹ N/m²

and that for water is 2.2 x 10⁹ N/m². The reciprocal of the bulk modulus is called the compressibility of the substance. The amount of compression of solids and liquids is seen to be very small.

The bulk modulus, B, is a material property that relates the change in volume with a change in pressure:

B = -V(¶P/¶V)_T = -(¶P/¶ln(V)) _T

A typical bulk modulus of medium density polyethylene is about 0.14 GPa (= 10⁹ N/m²). This is for a polyethylene that has a density of 0.93 g/cm³.

The bulk modulus can be calculated by distorting all of the dimensions of the unit cell using the program cellvec_xyz and calculating the energy as a function of the change in volume. This procedure uses the fact that the pressure is:

P = -(¶E/¶V)_T

and so

B = -V(¶²E/¶V²)_T

If we distort along all three dimensions of the unit cell by the same fractional amount we uniformly change the volume of the unit cell. If we plot the calculated energy as a function of the volume we will obtain a curve whose curvature (¶²E/¶V²) and unit cell volume V can be used to calculate the bulk modulus.

Young's modulus

For the description of the elastic properties of linear objects like wires, rods, columns, which are either stretched or compressed, a convenient parameter is the ratio of the stress to the strain, a parameter called the Young's modulus of the material. Young's modulus can be used to predict the elongation or compression of an object as long as the stress is less than the yield strength of the material.

Figure 1. Definition of terms used the calculation of elastic constants.

The relation between the external diagonal stress, P_i, and the diagonal components of the strain tensor, h_i, can be expressed as the following matrix equation:

To obtain the Young's modulus along the chain (assumed here to be in the x-direction according to our set-up of the unit cell above one needs to solve for h₁ at the condition P₁ = 1, P₂ = P₃ = 0 GPa.

The elastic constants are obtained from:

Where D_i and D_j are dimensionless displacements along unit cell vectors, V is the unit cell volume and U(D) and U(0) are the energies of the deformed and undeformed unit cell, respectively.

Example: Calculation of the Young's Modulus of Polyethylene

Using the program elastic_xyz.f we can create .car files with distortions along each of the unit directions, x, y and z. The program will prompt you for a four-letter root name and will then add "_x.car", "_y.car" and "_z.car", respectively, to create files:

root_x.car

root_y.car

root_z.car

with the same fractional distortion. We choose to sample the range from 0.96 to 1.04. We calculate the energy at each of the distorted unit cells using DFT methods. First, we perform a single point energy calculation. Plots are given for each of the dimensions.

Figure 2. The binding energy of the high density (orthorhombic) unit cell with displacements along the given unit cell directions for a single point energy calculation. The displacements are given in fractional coordinates along each unit cell direction.

Note that the minimum energy is not necessarily at the dimension given the x-ray crystal structure. A deviation likely arises from inaccuracy in the DFT calculation. The best model is obtained along the length of the carbon chains (x-direction). Here there is a clear minimum and Young's modulus could be calculated using these data. The deviation from experiment is only 1% along the x dimension.

Next we reoptimize the geometry for each deformed structure as shown in the Figure below.

Figure 3. The binding energy of the high density (orthorhombic) unit cell with displacements along the given unit cell directions for a reoptimized geometry at each deformation geometry. The displacements are given in fractional coordinates along each unit cell direction.

Along y and z there are clearly problems. The chain interactions along these dimensions are weak. DFT theory does not capture dispersion forces and so the weak chain interactions at larger distances may be an artifact.

Based on these calculated deformation energies we obtain the following deformation energies for D = 0.01 (using only one of the dimensions)

U(D_x) - U(0) = 1.0 ± 0.15, D = 0.01

U(D_y) - U(0) = 0.5 ± 0.25, D = 0.01

U(D_z) - U(0) = 0.25 ± 0.15 , D = 0.01

Ignoring off-diagonal terms in this simple example and substituting into the above formula (V = 2.534 x 4.93 x 7.4 = 92.4 ų) we find:

Y_x ~ 350 ± 54 , Y_x ~ 180 ± 90, Y_x ~ 90 ± 54 GPa

Experimental measurements of Y_x range from 230 to 360 depending on the method used. See reference 1.

There are at least two ways to improve on this result.

1. Sample smaller changes to determine whether there are well-defined minima. There could be more than one minimum given that there a number of interactions as the chains pushed together or pulled apart.

2. Better DFT calculations can be carried out by employing larger basis sets.

References

1. Miao et al. J. Chem. Phys. (2001), 115, 11317-11324