Concepts for Diagonalization of strain

This jupyter notebook is a worksheet that demostrates how coordinate transformation of strain works in the shear elastic modulus part works this program. You could clone this repo and change the key below and rerun it to see how the shear modulus solving mechanism is working for each of the elastic modulus

Finding the transformation matrix

Define a shear elastic modulus \(c_{ij}\) below as key

import numpy
from cij.util import c_

key = c_(46)

if not key.is_shear: raise RuntimeError(f"The strain c_{key} is not a shear modulus.")

Creating the corresponding fictitous strain \(e\) for \(c_{ij}\)

e = numpy.zeros((3, 3))

e[key.i[0] - 1, key.i[1] - 1] = 1
e[key.i[1] - 1, key.i[0] - 1] = 1

e[key.j[0] - 1, key.j[1] - 1] = 1
e[key.j[1] - 1, key.j[0] - 1] = 1

print(e)

[[0. 1. 0.]
 [1. 0. 1.]
 [0. 1. 0.]]

The transformation matrix \(T\) is the column-wise eigenvectors of the strain matrix \(e\)

evals, T = numpy.linalg.eig(e)
print(T)

[[ 5.00000000e-01  7.07106781e-01  5.00000000e-01]
 [-7.07106781e-01  9.02056208e-17  7.07106781e-01]
 [ 5.00000000e-01 -7.07106781e-01  5.00000000e-01]]

The \(e'\) in the rotated coordinate system is the matrix with eigenvalues of \(e\) in its diagonals

print(numpy.diag(evals))

[[-1.41421356e+00  0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  9.77950360e-17  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  1.41421356e+00]]

We can retrive the original strain with

\[e' = T^{-1} e T\]

print("e = ")
print(T @ numpy.diag(evals) @ T.T)

print("e = T e' T^{-1}:", numpy.allclose(e, T @ numpy.diag(evals) @ T.T))

e =
[[-1.11022302e-16  1.00000000e+00 -1.11022302e-16]
 [ 1.00000000e+00 -1.11022302e-16  1.00000000e+00]
 [-1.11022302e-16  1.00000000e+00  1.11022302e-16]]
e = T e' T^{-1}: True

or backwards

\[e = T e' T^{-1}\]

print("e' = ")
print(T.T @ e @ T)

print("e' = T^{-1} e T:", numpy.allclose(numpy.diag(evals), T.T @ e @ T))

e' =
[[-1.41421356e+00 -1.11022302e-16  0.00000000e+00]
 [-6.68036251e-17  4.93038066e-32  2.47214867e-16]
 [ 0.00000000e+00  3.33066907e-16  1.41421356e+00]]
e' = T^{-1} e T: True

Testing with symbols

import sympy

def round_expr(expr, num_digits=4):
    return expr.xreplace({n : round(n, num_digits) for n in expr.atoms(sympy.Number)})

The strain energy invariance equation

The strain energy \(E = \sum \frac{1}{2} c_{ijkl}e_{ij}e_{kl}\) should be same in the original and rotated coordinate system. Here we show the energy in both coordinate systems with the help of sympy. The function used to calculate the strain energy is

from cij.core.phonon_contribution import calculate_fictitious_strain_energy

The left hand side of the strain energy invariance equation, which is under the original coordinate system, the expression for the ficticious strain \(e\) defined above is

def resolve_elastic_modulus(key):
    return sympy.symbols("c_%d%d" % key.v)

calculate_fictitious_strain_energy(
    e,
    resolve_elastic_modulus
)

$\displaystyle 2.0 c_{44} + 4.0 c_{46} + 2.0 c_{66}$

and on the right hand side, which is under the rotated coordinate system, the expression for the rotated ficticious strain \(e'\) defined above is

def resolve_elastic_modulus_rotated(key):
    return sympy.symbols("c_%d'%d'" % key.v)

e_prime = T.T @ e @ T

calculate_fictitious_strain_energy(
    numpy.diag(numpy.diag(e_prime)),
    resolve_elastic_modulus_rotated,
)

$\displaystyle 1.0 c_{1'1'} - 2.0 c_{1'3'} + 1.0 c_{3'3'}$

The strain used for calculating strain-Grüneisen parameter from volume-Grüneisen parameter

a = sympy.symbols("e_11")
b = sympy.symbols("e_22")
c = sympy.symbols("e_33")

m = sympy.Matrix([
    [a, 0, 0],
    [0, b, 0],
    [0, 0, c]
])

round_expr(T.T @ m @ T)

$\displaystyle \left[\begin{matrix}0.25 e_{11} + 0.5 e_{22} + 0.25 e_{33} & 0.3536 e_{11} - 0.3536 e_{33} & 0.25 e_{11} - 0.5 e_{22} + 0.25 e_{33}\\0.3536 e_{11} - 0.3536 e_{33} & 0.5 e_{11} + 0.5 e_{33} & 0.3536 e_{11} - 0.3536 e_{33}\\0.25 e_{11} - 0.5 e_{22} + 0.25 e_{33} & 0.3536 e_{11} - 0.3536 e_{33} & 0.25 e_{11} + 0.5 e_{22} + 0.25 e_{33}\end{matrix}\right]$

a = sympy.symbols("α")
b = sympy.symbols("β")
c = sympy.symbols("γ")

m = sympy.Matrix([
    [1 + a, 0.   , 0.   ],
    [0.   , 1 + b, 0.   ],
    [0.   , 0.   , 1 + c]
])

round_expr(T.T @ m @ T)

$\displaystyle \left[\begin{matrix}0.25 α + 0.5 β + 0.25 γ + 1.0 & 0.3536 α - 0.3536 γ & 0.25 α - 0.5 β + 0.25 γ\\0.3536 α - 0.3536 γ & 0.5 α + 0.5 γ + 1.0 & 0.3536 α - 0.3536 γ\\0.25 α - 0.5 β + 0.25 γ & 0.3536 α - 0.3536 γ & 0.25 α + 0.5 β + 0.25 γ + 1.0\end{matrix}\right]$

Order for modulus calculation

from cij.util.voigt import _cij_sort_key
from IPython.display import display

sorted([
    c_(11), c_(12), c_(13), c_(22), c_(23), c_(33), c_(15), c_(25), c_(35), c_(46), c_(44), c_(55), c_(66)
], key=_cij_sort_key)

[11(1111),
 22(2222),
 33(3333),
 12(1122),
 13(1133),
 23(2233),
 44(2323),
 55(1313),
 66(1212),
 15(1113),
 25(2213),
 35(3313),
 46(2312)]