Post title: Cubic Elasticity

Post subtitle: Rotation of tensor of elastic constants

Created on: 18 Jun 2017

Updated on: 26 Dec 2024

Introduction

This post discusses the relationship between stresses \(\boldsymbol{\sigma}\) and strains \(\boldsymbol{\varepsilon}\) in crystals with cubic symmetry, specifically focusing on FCC and BCC crystal structures. Furthermore, the rotation of these tensors in three-dimensional space is illustrated.

The figure below shows a cubic crystal, which includes simple cubic (left), body-centered cubic (middle), and face-centered cubic (right) structures. Two coordinate systems are presented: the crystal coordinate system (\(x_c\), \(y_c\), \(z_c\)), aligned with the \(\langle100\rangle\) directions, and the global coordinate system (\(x\), \(y\), \(z\)), which is rotated relative to the crystal coordinate system about a vector \(\boldsymbol{o}\) by an angle \(\theta\).

The simplest constitutive model describing the relationship between stress (\(\boldsymbol{\sigma}\)) and strain (\(\boldsymbol{\varepsilon}\)) in these crystal structures is linear, commonly referred to as Hooke’s law. According to this law, stress is directly proportional to strain, therefore, we can write:

\[\begin{split} \sigma_{ij} &= \mathcal{D}_{ijkl} \varepsilon_{kl}\\ \boldsymbol{\sigma} &= \boldsymbol{D}\boldsymbol{\varepsilon} \end{split}, \label{}\tag{1}\]

where \(\boldsymbol{\sigma}\) is a second-rank stress tensor, \(\boldsymbol{\varepsilon}\) is a second-rank strain tensor, and \(\boldsymbol{D}\) is a fourth-rank tensor of elastic constants. This relationship can be written in the crystal coordinate system using the Voigt notation as follows:

\[\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}, \label{}\tag{2}\]

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 & D_{44}\end{bmatrix}\begin{bmatrix}\varepsilon_{xx} \\ \varepsilon_{yy} \\ \gamma_{xy}\end{bmatrix}. \label{}\tag{3}\]

It is important to note that the constitutive law is formulated using engineering strain. The shear strain components of the engineering strain, denoted as \(\gamma_{ij}\), are twice as large as the shear strain components of the tensorial strain, which can be expressed as \(\gamma_{ij}=2\varepsilon_{ij}\) for \(i \neq j\). Assuming small deformations, the tensorial strain can be calculated as \(\varepsilon_{ij}=1/2\left(\partial u_i/\partial x_j + \partial u_j/\partial x_i\right)\). Here, \(\vec{u}\) is the displacement field.

Depending on this assumption, if we consider plane strain conditions where \(\varepsilon_{zz}=0\), the constants of the strain tensor in two dimensions remain the same: \(C_{11}=D_{11}\), and \(C_{12}=D_{12}\). On the other hand, if we assume plane stress conditions, (\(\sigma_{zz}=0\)), \(C_{11}=\frac{D_{11}^2-D_{12}^2}{D_{11}}\), and \(C_{12}=\frac{D_{11}D_{12}-D_{12}^2}{D_{11}}\).

A useful parameter for describing cubic anisotropy is the Zener ratio, which is a dimensionless number calculated as follows: \(D_{44}/D'\) where \(D'=\frac{D_{11}-D_{12}}{2}\).

The elastic constants of selected materials, along with their corresponding Zener ratios, are summarized in the following table:

\(^a\) = (Simmons and Wang 1971), \(^b\) = (Jona and Marcus 2001)

Tensor rotation (rotation of elastic properties)

For interested readers, a more detailed analysis can be found in (Bao 2005).

This section begins by defining the rotation matrix We denote rotated tensors by a prime symbol (\('\)). A fourth-rank tensor is rotated as \(\boldsymbol{D}'=\boldsymbol{T}\cdot\boldsymbol{D}\cdot\boldsymbol{T}^T\), while a second-rank tensor, such as stress, is rotated as \(\boldsymbol{\sigma}' = \boldsymbol{R}\cdot\boldsymbol{\sigma}\cdot\boldsymbol{R}^T\), where \(\boldsymbol{R}\) is the rotation matrix.

Rotation (\(\boldsymbol{R}\)) and coordinate transformation (\(\boldsymbol{Q}\)) are distinct, with their relationship expressed as \(\boldsymbol{R} = \boldsymbol{Q}^T\). To begin, we consider the general form of a rotation matrix:

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

The \(\boldsymbol{T}\) matrix may then be assembled as

\[\boldsymbol{T} = \begin{bmatrix} \boldsymbol{T}^{(1)} & 2\boldsymbol{T}^{(2)} \\ \boldsymbol{T}^{(3)} & \boldsymbol{T}^{(4)} \\ \end{bmatrix}, \label{}\tag{5}\]

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}, } \label{}\tag{6}\]

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.. \label{}\tag{7}\]

Rotation about an axis \(\boldsymbol{o}=(o_1,o_2,o_3)\), \(|\boldsymbol{o}|=1\)

When rotation occurs about an arbitrary axis \(\boldsymbol{o}\), the rotation matrix \(\boldsymbol{R}\) simplifies to:

\[\boldsymbol{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} \label{}\tag{8}\]

Rotation of the Tensor of Elastic Constants \(\boldsymbol{C}\)

The constitutive law defines the tensor of elastic constants as a proportionality tensor relating stresses to engineering strains. The engineering strain shear terms \(\gamma_{ij}\) are twice as large as the tensorial strain shear terms \(\varepsilon_{ij}\). Considering this, the tensor of elastic constants is rotated using the following equation:

\[\boldsymbol{C}' = \boldsymbol{T}\boldsymbol{C}\boldsymbol{A}\boldsymbol{T}^{-1}\boldsymbol{A}^{-1} \label{}\tag{9}\]

where the matrix \(\boldsymbol{A}\) (also known as the Reuter’s matrix) is used to convert between the tensorial \(\boldsymbol{\varepsilon}\) and engineering \(\boldsymbol{\gamma}\) strains.

\[\boldsymbol{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} \label{}\tag{10}\]

The following Matlab script demonstrates how to rotate the tensor of elastic constants about an arbitrary vector \(\boldsymbol{o}\). This script uses 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 coordinate system (\(x\), \(y\), \(z\)).

% (c) 2018 Jakub Mikula
% PURPOSE: 
% Example of tensor rotation about a general axis
% Rotating tensor of elastic constants about an arbitrary vector o
%
% INPUT:
% 	> elastic constants D_11, D_12, D_44 
% 	> rotation axis o1, o2, o3 
% 	> rotation angle theta 
% OUTPUT:
%	> rotated tensor of elastic constants in matrix form D_rot

% -----------------------------------------------------------

% Material : [GPa]
	D_11 = 190
	D_12 = 161
%D_44 = (D_11-D_12)/2 % uncomment for isotropic elasticity
	D_44 = 42.3

% Axis to rotate about
% Make sure that |o|=1
	o1 = 0
	o2 = 0
	o3 = 1

% Specify the angle [RAD] 
	theta = -pi/4

% -----------------------------------------------------------

% Tensor of elastic constants (cubic elasticity) in the matrix form
% Voigth notation
D = [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];

% Matrix due to the conversion from engineering to tensorial strain	 
R(:,:) = 0.0d0;
R(1,1) = 1.0d0;
R(2,2) = 1.0d0;
R(3,3) = 1.0d0;
R(4,4) = 2.0d0;
R(5,5) = 2.0d0;
R(6,6) = 2.0d0;

invR(:,:) = R(:,:);
invR(4,4) = .5d0;
invR(5,5) = .5d0;
invR(6,6) = .5d0;


% Rotating about an arbitrary axis (between crystal and global coordinate system)
R_cg=[o1^2+(1-o1^2)*cos(theta) o1*o2*(1-cos(theta))-o3*sin(theta) o1*o3*(1-cos(theta))+o2*sin(theta)
   o1*o2*(1-cos(theta))+o3*sin(theta) o2^2+(1-o2^2)*cos(theta) o2*o3*(1-cos(theta))-o1*sin(theta)
   o1*o3*(1-cos(theta))-o2*sin(theta) o2*o3*(1-cos(theta))+o1*sin(theta) o3^2+(1-o3^2)*cos(theta)];


for i=1:3
for j=1:3

mmod = [1 1 1 1 1
		-1 2 2 2 2
		0 0 3 3 3
		1 1 1 4 4
		2 2 2 2 5];

 K1(i,j) = R_cg(i,j)^2;
 K2(i,j) = R_cg(i,mmod(j+1,3))*R_cg(i,mmod(j+2,3));
 K3(i,j) = R_cg(mmod(i+1,3),j)*R_cg(mmod(i+2,3),j);
 K4(i,j) = R_cg(mmod(i+1,3),mmod(j+1,3))*R_cg(mmod(i+2,3),mmod(j+2,3)) + ...
           R_cg(mmod(i+1,3),mmod(j+2,3))*R_cg(mmod(i+2,3),mmod(j+1,3));
end
end

T_cg(1:3,1:3) = K1;
T_cg(1:3,4:6) = 2.0d0*K2;
T_cg(4:6,1:3) = K3;
T_cg(4:6,4:6) = K4;

invT_cg(1:3,1:3) = K1';
invT_cg(1:3,4:6) = 2.0d0*K3';
invT_cg(4:6,1:3) = K2';
invT_cg(4:6,4:6) = K4';

% Rotated tensor of elastic constants
D_rot = T_cg*D*R*invT_cg*invR

Rotation of \(\boldsymbol{C}\) About \(\langle001\rangle\) in 2D

The rotation of the elastic constants can be derived using the following rotation matrices:

\[\boldsymbol{R}^{x_c} = \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos\theta & -\sin\theta \\ 0 & \sin\theta & \cos\theta \\ \end{bmatrix}, \boldsymbol{R}^{y_c} = \begin{bmatrix} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \\ \end{bmatrix}, \boldsymbol{R}^{z_C} = \begin{bmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}. \label{}\tag{11}\]

If we rotate the coordinate system about \(z_c\) (\([001]\)) by an angle \(\theta\), we can express the equations above in the following 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}. \label{}\tag{12}\]

The following Matlab code is provided:

% Purpose:  rotate tensor of elastic constants C into C_rot by angle theta
% Define:   theta [1]
%           C [3x3]
% Output:   C_rot [3x3]

C_rot(1,1) = C(1,1)*(cos(theta)^4+sin(theta)^4) + 2*C(1,2)*sin(theta)^2*cos(theta)^2 + C(3,3)*sin(2*theta)^2;
C_rot(1,2) = C(1,2)*(cos(theta)^4+sin(theta)^4) + 2*C(1,1)*sin(theta)^2*cos(theta)^2 - C(3,3)*sin(2*theta)^2;
C_rot(1,3) = 0.5*C(1,1)*sin(2*theta)*cos(2*theta) - 0.5*C(1,2)*sin(2*theta)*cos(2*theta) - 0.5*C(3,3)*sin(4*theta);
C_rot(2,1) = C_rot(1,2);
C_rot(2,2) = C_rot(1,1);
C_rot(2,3) = -C_rot(1,3);
C_rot(3,1) = C_rot(1,3);
C_rot(3,2) = -C_rot(1,3);
C_rot(3,3) = 0.5*C(1,1)*sin(2*theta)^2 - 0.5*C(1,2)*sin(2*theta)^2 + C(3,3)*cos(2*theta)^2;

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}. \label{}\tag{13}\]

The derivatives evaluated from the tensor of elastic constants vanish for isotropic materials, where \(C_{33}\) is related to \(C_{11}\) and \(C_{12}\) through the equation \(C_{33}=\frac{C_{11}-C_{12}}{2}\). This property has interesting implications for nanocrystalline materials (Mikula et al. 2019), where they are proportional to the bulk driving force contributing to grain boundary motion.

A Matlab code of the derivatives is provided below:

% Purpose:  rotate tensor of elastic constants C and calculate the derivative with respect to orientation theta into dC_rot
% Define:   theta [1]
%           C [3x3]
% Output:   dC_rot [3x3]

dC_rot(1,1) = -4*sin(theta)*cos(theta)*(2*cos(theta)^2*(C(1,1)-C(1,2)) - 4*cos(theta)^2*dC(3,3) - C(1,1) + C(1,2) + 2*C(3,3));
dC_rot(1,2) = 4*sin(theta)*cos(theta)*(2*cos(theta)^2*(C(1,1)-C(1,2)) - 4*cos(theta)^2*C(3,3) - C(1,1) + C(1,2) + 2*C(3,3));
dC_rot(1,3) = cos(2*theta)^2*(2*C(1,1)-2*C(1,2)) - 4*C(3,3)*cos(2*theta)^2 - C(1,1) + C(1,2) + 2*C(3,3);
dC_rot(2,1) = dC_rot(1,2);
dC_rot(2,2) = dC_rot(1,1);
dC_rot(2,3) = -dC_rot(1,3);
dC_rot(3,1) = dC_rot(1,3);
dC_rot(3,2) = -dC_rot(1,3);
dC_rot(3,3) = sin(4*theta)*(C(1,1) - C(1,2) - 2*C(3,3));

References