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?

matrix

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?

Stefano M · Accepted Answer

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 $6times 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 $6times 6$ matrices (containing products of directional cosines). Note that these matrices are not orthogonal; but invoking the invariance of $boldsymbol{sigma}^Tboldsymbol{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_sigmaboldsymbol{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.

Wolfgang Bangerth · Answer

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

Answered by Wolfgang Bangerth on July 31, 2020

user1318499 · Answer

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

debugging a rotation matrix for elastic constants

3 Answers

Recap theory

The answer

Things to be aware of

Add your own answers!

Ask a Question