Tensor rotation in Tilted Transverse Isotropy (TTI)

tti_anisotropy.py

"""
Tensor rotation in Tilted Transverse Isotropy (TTI) for the 36 stiffness elements
Discussion which I started here: https://www.researchgate.net/post/How_to_do_coordinate_transformation_of_a_6x6_stiffness_matrix_of_tilted_transverse_isotropic_medium

@author: Yohanes Nuwara
@email: ign.nuwara97@gmail.com

Possible extension of this work (may need collaboration):
* Verify that the rotation is right by calculating its tensor invariants
* Test in real lab: Does the theoretically rotated tensor represent the true stiffness of the TTI rock?
"""

A = 45.1
C = 34.03 
D = 8.28
M = 12.55 
F = 16.7
B = A - 2 * M

# stiffness tensor
S = np.array([[A, B, F, 0, 0, 0],
              [B, A, F, 0, 0, 0],
              [F, F, C, 0, 0, 0],
              [0, 0, 0, D, 0, 0],
              [0, 0, 0, 0, D, 0],
              [0, 0, 0, 0, 0, M]])

" Strike and dip data "

strike = np.arange(0, 180, 5) # azimuth from north
dip = 30
theta = strike 
psi = dip 

c11 = []; c12 = []; c13 = []
c21 = []; c22 = []; c23 = []
c31 = []; c32 = []; c33 = []

c14 = []; c15 = []; c16 = []
c24 = []; c25 = []; c26 = []
c34 = []; c35 = []; c36 = []

c41 = []; c42 = []; c43 = []
c51 = []; c52 = []; c53 = []
c61 = []; c62 = []; c63 = []

c44 = []; c45 = []; c46 = []
c54 = []; c55 = []; c56 = []
c64 = []; c65 = []; c66 = []

for i in range(len(strike)):
  " direction cosines of rotation along Z axis (rotate strike) "

  theta = strike[i] # transformation angle

  a1 = np.cos(np.deg2rad(theta)); a2 = -np.sin(np.deg2rad(theta)); a3 = 0
  b1 = np.sin(np.deg2rad(theta)); b2 = np.cos(np.deg2rad(theta)); b3 = 0
  c1 = 0; c2 = 0; c3 = 1

  Z = np.array([[a1, a2, a3],
                [b1, b2, b3],
                [c1, c2, c3]])
  
  " direction cosines of rotation along X axis (rotate dip) "

  psi = dip # transformation angle

  d1 = 1; d2 = 0; d3 = 0
  e1 = 0; e2 = np.cos(np.deg2rad(psi)); e3 = -np.sin(np.deg2rad(psi))
  f1 = 0; f2 = np.sin(np.deg2rad(psi)); f3 = np.cos(np.deg2rad(psi))

  X = np.array([[d1, d2, d3],
                [e1, e2, e3],
                [f1, f2, f3]])

  " direction cosine of rotation combination of strike and dip "

  L = np.dot(X, Z)

  " direction cosine elements "

  l1 = L[0][0]; l2 = L[0][1]; l3 = L[0][2]
  m1 = L[1][0]; m2 = L[1][1]; m3 = L[1][2]
  n1 = L[2][0]; n2 = L[2][1]; n3 = L[2][2]

  " rotation tensor "

  Y = np.array([[l1**2,   m1**2,   n1**2,   2*m1*n1,        2*n1*l1,        2*l1*m1],
                [l2**2,   m2**2,   n2**2,   2*m2*n2,        2*n2*l2,        2*l2*m2],
                [l3**2,   m3**2,   n3**2,   2*m3*n3,        2*n3*l3,        2*l3*m3],
                [l2*l3,   m2*m3,   n2*n3,   m2*n3 + m3*n2,  n2*l3 + n3*l2,  m2*l3 + m3*l2],
                [l3*l1,   m3*m1,   n3*n1,   m3*n1 + m1*n3,  n3*l1 + n1*l3,  m3*l1 + m1*l3],
                [l1*l2,   m1*m2,   n1*n2,   m1*n2 + m2*n1,  n1*l2 + n2*l1,  m1*l2 + m2*l1]])

  " transformation of stiffness "

  Y_inv = np.linalg.inv(Y)
  S_dot = np.dot(Y, S)
  S_trans = np.dot(S_dot, Y_inv)

  " results of the stiffness elements "

  c11.append(float(S_trans[0][0])); c12.append(float(S_trans[0][1])); c13.append(float(S_trans[0][2]))
  c21.append(float(S_trans[1][0])); c22.append(float(S_trans[1][1])); c23.append(float(S_trans[1][2]))
  c31.append(float(S_trans[2][0])); c32.append(float(S_trans[2][1])); c33.append(float(S_trans[2][2]))

  c14.append(float(S_trans[0][3])); c15.append(float(S_trans[0][4])); c16.append(float(S_trans[0][5]))
  c24.append(float(S_trans[1][3])); c25.append(float(S_trans[1][4])); c26.append(float(S_trans[1][5]))
  c34.append(float(S_trans[2][3])); c35.append(float(S_trans[2][4])); c36.append(float(S_trans[2][5]))

  c41.append(float(S_trans[3][0])); c42.append(float(S_trans[3][1])); c43.append(float(S_trans[3][2]))
  c51.append(float(S_trans[4][0])); c52.append(float(S_trans[4][1])); c53.append(float(S_trans[4][2]))
  c61.append(float(S_trans[5][0])); c62.append(float(S_trans[5][1])); c63.append(float(S_trans[5][2]))

  c44.append(float(S_trans[3][3])); c45.append(float(S_trans[3][4])); c46.append(float(S_trans[3][5]))
  c54.append(float(S_trans[4][3])); c55.append(float(S_trans[4][4])); c56.append(float(S_trans[4][5]))
  c64.append(float(S_trans[5][3])); c65.append(float(S_trans[5][4])); c66.append(float(S_trans[5][5]))

plt.plot(strike, c22, '.-')
plt.title('Anisotropy Stiffness Element C22', size=17, pad=10)
plt.xlabel('Strike of bedding plane ($N...^oE$)', size=15)
plt.ylabel('Stiffness (GPa)', size=15)
plt.xlim(0, max(strike))

One of the elements (C22) stiffness after rotation: