debugging a rotation matrix for elastic constants

I have a transformation matrix which takes in the elastic constants from the local $rtl$ coordinates and then converts the elastic constants to the global $xyz$ coordinates via a rotation about the $z$ axis (i.e., $l$ doesn't change with respect to $z$, at lest at this stage). From what I gather the standard way to do this is to create a $6 \times 6$ matrix of directional cosines, along with a $6 \times 6$ matrix of elastic constants (the stiffness matrix) and then use $G^T C G$. where $G$ is the rotation matrix and $C$ is the stiffness matrix. I did this, and then to validate it I used the analytical solutions in Lekhnitskii's text 'Theory of elasticity of an anisotropic elastic body' for a cylinder rotated around the $z$ axis. To check I had coded Lekhnitskii's matrix correctly I used the 5 rotational invariants presented partially in his book and elsewhere on the web. Now his matrix is working (at lest as far as the 5 invariants hold), however the $G^T C G$ system does not match with his solution, (and most of the invariants don't hold for my $G^T C G$ system either).

I have two questions: Is there a more complete way to test that Lekhnitskii's matrix is working correctly?

And is there a systematic way to debug my $G^T C G$ system to find out why it is not working?

• I would like to follow-up on previous answers. I understand that the expressions for T are terrible. Is there a good and reliable source where I can look these expressions up? – Sagar Shah Jul 21 '20 at 12:49
• @SagarShah, if you are interested you should ask a new question. – nicoguaro Jul 21 '20 at 14:28

You are facing one of the most tedious and error prone aspect of elasticity theory: change of reference frame in engineering (or Voigt) notation.

Recap theory

If $\boldsymbol{\epsilon}$ and $\boldsymbol{\sigma}$ are the strain and stress engineering components (represented as $6\times 1$ vectors), then you can change reference of frame with the following transformation: \begin{gather} \boldsymbol{\epsilon}' = T_\epsilon \, \boldsymbol{\epsilon}\\ \boldsymbol{\sigma}' = T_\sigma \, \boldsymbol{\sigma} \end{gather} where $T_\sigma, T_\epsilon$ are $6\times 6$ matrices (containing products of directional cosines). Note that these matrices are not orthogonal; but invoking the invariance of $\boldsymbol{\sigma}^T\boldsymbol{\epsilon}$ we have $$T_\epsilon = T_\sigma ^ {-T}$$ For the stifness matrix $D$, remembering that $\boldsymbol{\sigma} = D \boldsymbol{\epsilon}$ becomes $\boldsymbol{\sigma}' = D' \boldsymbol{\epsilon}'$ we have $$D' = T_\sigma \, D\, T_\epsilon^{-1} = T_\sigma \, D\, T_\sigma^T$$ Note however that the compliance matrix $C = D^{-1}$ follows a different rule: $$C\,{}' = T_\epsilon \, C\, T_\sigma^{-1} = T_\epsilon \, C\, T_\epsilon^T$$ Expressions for $T_\epsilon$ and $T_\sigma$ are terrible and forgive me if I do not copy them here.

The easiest way to check for a correct expression of $T_\sigma$ is to compute the components of $\boldsymbol{\sigma}'$ in tensor notation $\sigma_{ij}'=a_{im}a_{jl}\sigma_{ml}$ where $a_{ij}$ are the directional cosines of the given transformation. In matrix notation this amounts to $\Sigma' = A \Sigma A^T$, where $A$ is the usual unitary rotation matrix, so its very hard to do a mistake here. Then you compare the components of $T_\sigma\boldsymbol{\sigma}$ and $\Sigma'$ to see if they match. Test this for a number of arbitrary rotations.

Things to be aware of

The expression of $T_\sigma$ and $T_\epsilon$ depend on the ordering of $\boldsymbol{\sigma}$ and $\boldsymbol{\epsilon}$: double check that you are consistent.

Compliance and stiffness matrix in engineering notation have different transformation laws, owing to the fact that $T$ matrices are not unitary.

$T_\sigma$ and $T_\epsilon$ are almost equal, up to a factor of two, because engineering shear strains are twice the corresponding tensor components.

Here's a clearer way to transform a 6x6 stiffness matrix - treat it as a 3x3x3x3 tensor:

   For k = 1 To 3
For l = 1 To 3
For s = 1 To 3
For t = 1 To 3

rotatedTensor(k,l,s,t) = 0

For m = 1 To 3
For n = 1 To 3
For p = 1 To 3
For r = 1 To 3

rotatedTensor(k,l,s,t) += tensor(m,n,p,r) * T(k, m) * T(l, n) * T(s, p) * T(t, r)

Next
Next
Next
Next

Next
Next
Next
Next


T is an ordinary 3x3 rotation matrix. The input is tensor and the output is rotatedtensor. Wherever 4 indices appear, convert them to the 2-index form used in the stiffness matrix. That way you can store the input and output as 6x6 matrices and just use the 4 indices to make the code more readable.

The pseudocode above is independent of the particular mapping from 4-indices to 2-indices. You have to apply that mapping yourself at each point where it refers to a tensor using 4 indices. For example, you might define a function:

double tensor(a,b,c,d){
if(a==1 and b==1 and c==1 and d==1){
return matrix(1,1)
}
if(a==1 and b==1 and c==2 and d==2){
return matrix(1,2)
}
if(a==1 and b==1 and c==3 and d==3){
return matrix(1,3)
}
if(a==1 and b==1 and c==2 and d==3){
return matrix(1,4)
}
etc ...
}


As a matter of general truth: Life becomes a lot easier and less error prone if you use $3\times 3$ matrices instead of writing tensors as 6-component vectors. The laws of transforming $3\times 3$ are much more obvious than the transformation of 6-component vector representation of these symmetric matrices.

• Unfortunately most of the engineering world uses Voigt notation... In the specific case of anisotropic material behavior it is very handy to have the forth order symmetric stiffness tensor $D_{ijhk}$ represented as a $6\times 6$ symmetric matrix $D$. And the OP question was on change of reference of the stiffness tensor/matrix. – Stefano M Jul 31 '13 at 16:40
• I know that engineers like to use this notation (despite the fact that I think that it's a lot more difficult to use). All I wanted to point out is that the question becomes simple if you wrote it in tensor form. One way to make the problem simple is if you know how to transform back and forth between the two notations. – Wolfgang Bangerth Jul 31 '13 at 19:59