TransWikia.com

debugging a rotation matrix for elastic constants

Computational Science Asked on July 31, 2020

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?

3 Answers

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.

Correct answer by Stefano M on July 31, 2020

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

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

Answered by user1318499 on July 31, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP