What is cubic symmetry? How to rotate the corresponding tensor of elastic constants?
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$.
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,
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)
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):
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} $$
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
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,
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$.
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