Let's solve the Christoffel equation system using Sympy.
First, let's build a cij matrix using Sympy.
import sympy
def make_cij_matrix_isotropic():
c11, c44 = sympy.symbols('c_11 c_44')
return sympy.Matrix([[c11, c11- 2*c44, c11- 2*c44, 0, 0, 0],
[c11- 2*c44, c11, c11- 2*c44, 0, 0, 0],
[c11- 2*c44, c11- 2*c44, c11, 0, 0, 0],
[0, 0, 0, c44, 0, 0],
[0, 0, 0, 0 , c44, 0],
[0, 0, 0, 0, 0, c44]])
cij = make_cij_matrix_isotropic()
cij
def to_cijkl(cij):
"""Converts the Voigt-notation cij to a 4th order tensor cijkl."""
coord_mapping = {(1, 1): 1,
(2, 2): 2,
(3, 3): 3,
(2, 3): 4,
(1, 3): 5,
(1, 2): 6,
(2, 1): 6,
(3, 1): 5,
(3, 2): 4}
cijkl = sympy.MutableDenseNDimArray([0] * 81, shape=(3, 3, 3, 3))
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
u = coord_mapping[(i + 1, j + 1)]
v = coord_mapping[(k + 1, l + 1)]
cijkl[i, j, k, l] = cij[u - 1, v - 1]
return cijkl
cijkl = to_cijkl(cij)
cijkl
rho = sympy.symbols('rho')
def assemble_christoffel_matrix(cijkl, q, rho):
M = sympy.zeros(3, 3)
for i in range(3):
for j in range(3):
for n in range(3):
for m in range(3):
M[i, j] += q[n] * cijkl[i, n, m, j] * q[m]
return M / rho
q = sympy.Matrix([0, 0, 1]).reshape(3, 1)
assemble_christoffel_matrix(cijkl, q, rho)
q = sympy.Matrix([1, 0, 0]).reshape(3, 1)
assemble_christoffel_matrix(cijkl, q, rho)
q = sympy.Matrix([0, 1, 0]).reshape(3, 1)
assemble_christoffel_matrix(cijkl, q, rho)
Let's make a vector of squared speeds:
vs_1, vs_2, vs_3 = sympy.symbols('vs_1 vs_2 vs_3')
squared_speeds = sympy.Matrix([vs_1, vs_2, vs_3]).reshape(3, 1)
squared_speeds
q = sympy.Matrix([0, 1, 0]).reshape(3, 1)
M = assemble_christoffel_matrix(cijkl, q, rho)
M - sympy.diag(*[vs_1, vs_2, vs_3])
eigenvects = (M - sympy.diag(*[vs_1, vs_2, vs_3])).eigenvects()
eigenvects
[(c_11/rho - vs_2, 1, [Matrix([ [0], [1], [0]])]), (c_44/rho - vs_1, 1, [Matrix([ [1], [0], [0]])]), (c_44/rho - vs_3, 1, [Matrix([ [0], [0], [1]])])]
Let's print some results using HTML:
from IPython.display import HTML, display
def pprint_eigenvects(wave_vector, eigenvects):
html = ''
html += f'<b>incidence : {wave_vector._repr_latex_()}</b>'
for eigenval, mult, vect in eigenvects:
for var in [vs_1, vs_2, vs_3]:
if len(eigenval.find(var)) > 0:
sol = sympy.solve(sympy.Eq(eigenval, 0), var)
html += f'vitesse {sympy.sqrt(sol[0])._repr_latex_()}'
html += f' selon la polarisation {vect[0]._repr_latex_()}, '
html += '<br>'
return HTML(html)
for direction in range(3):
wave_vector = sympy.Matrix([0, 0, 0])
wave_vector[direction] = 1
M = assemble_christoffel_matrix(cijkl, wave_vector, rho)
eigenvects = (M - sympy.diag(*[vs_1, vs_2, vs_3])).eigenvects()
display(pprint_eigenvects(wave_vector, eigenvects))
def make_cij_matrix_cubic():
c11, c12, c44 = sympy.symbols('c_11 c_12 c_44')
return sympy.Matrix( [[c11, c12, c12, 0, 0, 0],
[c12, c11, c12, 0, 0, 0],
[c12, c12, c11, 0, 0, 0],
[0, 0, 0, c44, 0, 0],
[0, 0, 0, 0 , c44, 0],
[0, 0, 0, 0, 0, c44]])
cij = make_cij_matrix_cubic()
cij
cijkl = to_cijkl(cij)
for direction in range(3):
wave_vector = sympy.Matrix([0, 0, 0])
wave_vector[direction] = 1
M = assemble_christoffel_matrix(cijkl, wave_vector, rho)
eigenvects = (M - sympy.diag(*[vs_1, vs_2, vs_3])).eigenvects()
display(pprint_eigenvects(wave_vector, eigenvects))
def make_cij_matrix_orthotropic():
c11, c12, c13, c22, c23, c33, c44, c55, c66 = sympy.symbols('c_11 c_12 c_13 c_22 c_23 c_33 c_44 c_55 c_66')
return sympy.Matrix( [[c11, c12, c13, 0, 0, 0],
[c12, c22, c23, 0, 0, 0],
[c13, c23, c33, 0, 0, 0],
[0, 0, 0, c44, 0, 0],
[0, 0, 0, 0 , c55, 0],
[0, 0, 0, 0, 0, c66]])
cij = make_cij_matrix_orthotropic()
cij
cijkl = to_cijkl(cij)
for direction in range(3):
wave_vector = sympy.Matrix([0, 0, 0])
wave_vector[direction] = 1
M = assemble_christoffel_matrix(cijkl, wave_vector, rho)
eigenvects = (M - sympy.diag(*[vs_1, vs_2, vs_3])).eigenvects()
display(pprint_eigenvects(wave_vector, eigenvects))