debugging a rotation matrix for elastic constants

Question

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?

Answer 1

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 \begin{equation} T_\epsilon = T_\sigma ^ {-T} \end{equation} For the stifness matrix $D$, remembering that $\boldsymbol{\sigma} = D \boldsymbol{\epsilon}$ becomes $\boldsymbol{\sigma}' = D' \boldsymbol{\epsilon}'$ we have \begin{equation} D' = T_\sigma \, D\, T_\epsilon^{-1} = T_\sigma \, D\, T_\sigma^T \end{equation} Note however that the compliance matrix $C = D^{-1}$ follows a different rule: \begin{equation} C\,{}' = T_\epsilon \, C\, T_\sigma^{-1} = T_\epsilon \, C\, T_\epsilon^T \end{equation} Expressions for $T_\epsilon$ and $T_\sigma$ are terrible and forgive me if I do not copy them here.

The answer

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.

Answer 2

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 ...
}

Answer 3

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.