Cubic Elasticity

What is cubic symmetry? How to rotate the corresponding tensor of elastic constants?

Introduction

In this post we establish a simple relationship between stresses $\boldsymbol{\sigma}$ and strains $\boldsymbol{\varepsilon}$ of crystals with cubic symmetry (FCC, BCC). Furthermore, we illustrate rotation of the tensors.

The following figure demonstrates a cubic crystal (simple cubic, body-centered cubic, and face-centered cubic) with two coordinate systems: crystal coordinate system ($x_c$,$y_c$,$z_c$) aligned with the <100> directions and a general global coordinate system ($x$, $y$, $z$) chosen by convenience and rotated from the crystal coordinate system about a vector $\boldsymbol{o}$ by angle $\theta$.

Fig.1 Example of a single crystal: crystal ($x'$,$y'$,$z'$) and global ($\mathbf{x}$,$\mathbf{y}$,$\mathbf{z}$) coordinate system.

The simplest constitutive model (relation between stress $\mathbf{\sigma}$ and strain $\mathbf{\varepsilon}$) is linear. Also known as the Hook’s law, the stress is linearly proportional to strain:

\[\boldsymbol{\sigma}=\boldsymbol{D}\boldsymbol{\varepsilon},\]

where $\boldsymbol{\sigma}$ is a second-rand stress tensor, $\boldsymbol{\varepsilon}$ is a second-rank strain tensor, and $\boldsymbol{D}$ is a fourth-rank tensor of elastic constants. For a three-dimensional stress state in cartesian coordinates and using the Voight notation, this relation can be in the crystal coordinate system written as

$$\begin{bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{zz} \\ \sigma_{yz} \\ \sigma_{zx} \\ \sigma_{xy} \end{bmatrix} = \begin{bmatrix} D_{11} & D_{12} & D_{12} & 0 & 0 & 0 \\ D_{12} & D_{11} & D_{12} & 0 & 0 & 0 \\ D_{12} & D_{12} & D_{11} & 0 & 0 & 0 \\ 0 & 0 & 0 & D_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & D_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & D_{44} \end{bmatrix} \begin{bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ \varepsilon_{zz} \\ \gamma_{yz} \\ \gamma_{zx} \\ \gamma_{xy} \end{bmatrix},$$

which in 2D simplifies into

\[\begin{bmatrix}\sigma_{xx} \\ \sigma_{yy} \\ \tau_{xy}\end{bmatrix} = \begin{bmatrix}C_{11} & C_{12} & 0 \\ C_{12} & C_{11} & 0 \\ 0 & 0 & C_{33}\end{bmatrix}\begin{bmatrix}\varepsilon_{xx} \\ \varepsilon_{yy} \\ \gamma_{xy}\end{bmatrix}.\]

The constitutive law is given considering the engineering strain. The shear strain components of the engineering strain $\gamma_{ij}$ are twice that much of the shear strain components of the tensorial strain $\gamma_{ij}=2\varepsilon_{ij}$. Assumming small deformations, the tensorial strain is calculated as $\varepsilon_{ij}=1/2\left(\partial u_i/\partial x_j + \partial u_j/\partial x_i\right)$.

Now, depending on this assumption, if we consider plane strain conditions ($\varepsilon_{zz}=0$), the constants of the strain tensor in 2D remain the same $C_{11}=D_{11}$ , $C_{12}=D_{12}$, and $C_{33}=D_{44}$, otherwise if we assume plane stress conditions ($\sigma_{zz}=0$), $C_{11}=\frac{D_{11}^2-D_{12}^2}{D_{11}}$, $C_{12}=\frac{D_{11}D_{12}-D_{12}^2}{D_{11}}$, $C_{33}=D_{44}$.

A useful parameter to describe the cubic anisotropy is using the Zener ratio. This non-dimensional number is calculated as $D_{44}/D’$ where $D’=\frac{D_{11}-D_{12}}{2}$.

The following table summarizes the elastic constants of selected materials together with the Zener ratio,

Table 1: Elastic constants of some FCC metals
Metal $C_{11}$ (GPa) $C_{12}$ (GPa) $C_{44}$ (GPa) $C_{44}/C’$
Pb$^a$ 55.6 45.4 19.4 3.8
Ag$^a$ 131.5 97.3 51.1 2.99
Au$^a$ 201 170 46 2.97
Cu$^b$ 225 153 115 3.19
Ni$^a$ 261 151 132 2.4
Al$^a$ 114 62 32 1.23
Pd$^a$ 232 176 71.2 2.45
Pt$^a$ 358 253 77.5 1.47

$^a$ Simmons G, Wang H. Single crystal elastic constants and calculated aggregate properties: a handbook. Cambridge, MA: MIT; 1971

$^b$ Jona F, Marcus PM. Phys Rev B 2001;63:094113

The table has been adapted from: S. Kibey et al. Predicting twinning stress in fcc metals: Linking twin-energy pathways to twin nucleation, Acta Materialia 55 (2007)



Tensor Rotation (Rotation of Elastic Properties)

In this section, we denote rotated tensors by $’$. A fourth-rank tensor is rotated as $\mathbf{D}’=\mathbf{T}\cdot\mathbf{D}\cdot\mathbf{T}^T$ and a second-rank tensor, such as stress, as $\mathbf{\sigma}’ = \mathbf{R}\cdot\mathbf{\sigma}\cdot\mathbf{R}^T$, where $\mathbf{R}$ is the rotation matrix.

We distinguish between rotation $\mathbf{R}$ and coordinate transformation $\mathbf{Q}$. It holds that $\mathbf{R}=\mathbf{Q}^T$. We begin with the general form of rotation matrix

\[\mathbf{R} = \begin{bmatrix} \cos(\boldsymbol{x}',\boldsymbol{x}) & \cos(\mathbf{y}',\mathbf{x}) & \cos(\mathbf{z}',\mathbf{x}) \\ \cos(\boldsymbol{x}',\boldsymbol{y}) & \cos(\mathbf{y}',\mathbf{y}) & \cos(\mathbf{z}',\mathbf{y}) \\ \cos(\boldsymbol{x}',\boldsymbol{z}) & \cos(\mathbf{y}',\mathbf{z}) & \cos(\mathbf{z}',\mathbf{z}) \\ \end{bmatrix}.\]

The $\mathbf{T}$ matrix may then be assembled as

\[\mathbf{T} = \begin{bmatrix} \mathbf{T}^{(1)} & 2\mathbf{T}^{(2)} \\ \mathbf{T}^{(3)} & \mathbf{T}^{(4)} \\ \end{bmatrix},\]

where

$$ {\small \begin{split} T_{ij}^{(1)} &= R_{ij}^2 \\ T_{ij}^{(2)} &= R_{i \textrm{mod}(j+1,3)}R_{i \textrm{mod}(j+2,3)} \\ T_{ij}^{(3)} &= R_{\textrm{mod}(i+1,3)j}R_{\textrm{mod}(i+2,3)j} \\ T_{ij}^{(4)} & = R_{\textrm{mod}(i+1,3)\textrm{mod}(j+1,3)}R_{\textrm{mod}(i+2,3)\textrm{mod}(j+2,3)} + R_{\textrm{mod}(i+1,3)\textrm{mod}(j+2,3)}R_{\textrm{mod}(i+2,3)\textrm{mod}(j+1,3)} \\ \end{split}, } $$

where the modulo function is defined as

\[\textrm{mod(i,3)} = \left\{ \begin{array}{c c} i & i\leq 3 \\ i-3 & i>3 \\ \end{array} \right..\]

This can be easily programmed (a Matlab code example provided below):

Rotation About an Axis $\mathbf{o}=(o_1,o_2,o_3)$, $|\mathbf{o}|=1$

In a special case, when rotating about an arbitrary axis $\mathbf{o}$, the rotation matrix $\mathbf{R}$ can be simplified into

$$ \mathbf{R} = \begin{bmatrix} o_1^2 + (1-o1^2)\cos\theta & o_1o_2(1-\cos\theta)-o_3\sin\theta & o_1o_3(1-\cos\theta)+o_2\sin\theta \\ o_1o_2(1-\cos\theta)+o_3\sin\theta & o_2^2+(1-o_2^2)\cos\theta & o_2o_3(1-\cos\theta)-o_1\sin\theta \\ o_1o_3(1-\cos\theta)-o_2\sin\theta & o_2o_3(1-\cos\theta) + o_1\sin\theta & o_3^2+(1-o_3^2)\cos\theta \\ \end{bmatrix} $$

Rotation the Tensor of Elastic Constants $\mathbf{C}$

Remember that the constitutive law is writen such that the tensor of elastic constants is a proportionality tensor between stresses and engineering strains. The engineering strain shear terms $\gamma_{ij}$ are twice that much as the tensorial strain shear terms $\varepsilon_{ij}$. Taking this into account, the tensor of elastic constants is rotated as

\[\mathbf{C}' = \mathbf{T}\mathbf{C}\mathbf{A}\mathbf{T}^{-1}\mathbf{A}^{-1}\]

where the matrix $\mathbf{A}$ converts between the rensorial $\mathbf{\varepsilon}$ and engineering $\mathbf{\gamma}$ strains

\[\mathbf{A} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 2 & 0 & 0\\ 0 & 0 & 0 & 0 & 2 & 0\\ 0 & 0 & 0 & 0 & 0 & 2\\ \end{bmatrix}\]

Rotating of the tensor of elastic constants about an arbitrary vector $\mathbf{o}$ is provided in the Matlab script below. This script demonstrates the background image of this post, in which the elastic constants $D_{11}$, $D_{12}$, and $D_{44}$ are given in the crystal coordinate system, and the tensor of elastic constants is rotated into the global coodinate system ($x$, $y$, $z$).

Interactive form available at: http://jakubmikula.com/interactives/rotate_rank4.html

Rotation of $\mathbf{C}$ About <001> in 2D

The rotation of the elastic elements can be derived using the following rotation matrices,

\(\mathbf{R}^{x_c} = \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos\theta & -\sin\theta \\ 0 & \sin\theta & \cos\theta \\ \end{bmatrix},\) \(\mathbf{R}^{y_c} = \begin{bmatrix} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \\ \end{bmatrix},\) \(\mathbf{R}^{z_C} = \begin{bmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}.\)

Rotating about $z_c$ ($[001]$) by an angle $\theta$, the equations above can be then written in the form

\[\begin{array}{l} C_{11}' = C_{11}\left(\cos^4\theta + \sin^4\theta\right) + 2C_{12}\sin^2\theta\cos^2\theta + C_{33}\sin^22\theta \\ C_{12}' = C_{12}\left(\cos^4\theta + \sin^4\theta\right) + 2C_{11}\sin^2\theta\cos^2\theta - C_{33}\sin^22\theta \\ C_{13}' = \frac{1}{2}C_{11}\sin 2\theta \cos 2\theta - \frac{1}{2}C_{12}\sin 2\theta \cos 2\theta - \frac{1}{2}C_{33}\sin 4\theta \\ C_{21}' = C_{12}' \\ C_{22}' = C_{11}' \\ C_{23}' = -C_{13}' \\ C_{31}' = C_{13}' \\ C_{32}' = -C_{13}' \\ C_{33}' = \frac{1}{2}C_{11}\sin^2 2\theta - \frac{1}{2}C_{12} \sin^2 2\theta + C_{33}\cos^2 2\theta \\ \end{array}.\]

Below is provided a Matlab code,

The derivatives of these with respect to the orientation $\theta$ are,

$$ \begin{array}{l} \frac{\partial C_{11}'}{\partial \theta} = - 4\sin\theta\cos\theta \left(2\cos^2\theta(C_{11}-C_{12})-4\cos^2\theta C_{33} - C_{11} + C_{12} + 2C_{33}\right)\\ \frac{\partial C_{12}'}{\partial \theta} = 4\sin\theta\cos\theta\left(2\cos^2\theta(C_{11}-C_{12})-4\cos^2\theta C_{33} - C_{11} + C_{12} + 2C_{33}\right)\\ \frac{\partial C_{13}'}{\partial \theta} = \cos^22\theta(2C_{11}-2C_{12})-4C_{33}\cos^22\theta - C_{11} + C_{12} + 2C_{33}\\ \frac{\partial C_{33}'}{\partial \theta} = \sin4\theta(C_{11}-C_{12}-2C_{33}) \end{array}. $$

The derivatives yield zero for isotropic elasticity ($C_{33}=\frac{C_{11}-C_{12}}{2}$). In my recent publication (Mikula et al.,SMS,2018) I discuss how these are proportional to the driving force to move grain boundaries in nanocrystalline materials.

Below is provided a Matlab code,



Stress-Strain Curves

As a material model example we choose gold and aluminium nanocrystals. Gold is relatively anisotropic compared to aluminium. The stress-strain curves $\sigma_{xx}$ vs. $\varepsilon_{xx}$ of the two crystals loaded under different orientations can be obtained as a line with the slope equal to the ratio $\sigma_{xx}/\varepsilon_{xx}$. The following figure demonstrates this ratio (stiffness) for different crystal orientations $\theta$ when rotating the single crystal of aluminium and gold about $z_c$.

Fig.2 Ratio $\sigma_{xx}/\varepsilon_{xx}$ for single crystals of Au and Al loaded in different directions rotating about $z_c$.

Sources

Bao, M. (2005). Analysis and design principles of MEMS devices. Amsterdam;Oxford;: Elsevier. (in Chapter 6).

http://www.mse.mtu.edu/~drjohn/my4150/compositesdesign/cd3/cd1.html

http://solidmechanics.org/text/Chapter3_2/Chapter3_2.htm

https://www.brown.edu/Departments/Engineering/Courses/EN224/anis_general/anis_general.htm