# SymmetrizedArray of stiffness/compliance tensor

The stiffness and compliance tensors of a material are 3x3x3x3 tensors relating stress and strain.

 (* a general 3x3x3x3 tensor *)
general = SymmetrizedArray[
pos_ :> Subscript[c, StringJoin[ToString /@ pos]],
{3, 3, 3, 3},
{}];


Since stress and strain are each symmetric, and the energy is a quadratic form, we can reduce the number of independent components to 21 for an anisotropic material.

(* an anisotropic stiffness tensor *)
anisotropic = SymmetrizedArray[
pos_ :> Subscript[c, StringJoin[ToString /@ pos]],
{3, 3, 3, 3},
{(* stress (strain) is symmetric *){Cycles[{{1, 2}}], 1},
(* strain (stress) is symmetric *){Cycles[{{3, 4}}], 1},
(* energy is quadratic *) {Cycles[{{1, 3}, {2, 4}}], 1}
}];


Material symmetries (invariance under various rotations and inversions) can further reduce the number of independent components.

How can we tell Mathematica about the material symmetries?

(* an orthotropic stiffness tensor *)
orthotropic = SymmetrizedArray[
pos_ :> Subscript[c, StringJoin[ToString /@ pos]],
{3, 3, 3, 3},
{{Cycles[{{1, 2}}], 1},
{Cycles[{{3, 4}}], 1},
{Cycles[{{1, 3}, {2, 4}}], 1},
(* WHAT GOES HERE? *)
}];


I tried using TensorSymmetry on an orthotropic tensor but it just returned the anisotropic cycles. The solution should give the same result as the following:

orthotropicRules = {
Subscript[c, "1123"] -> 0, Subscript[c, "1113"] -> 0, Subscript[c, "1112"] -> 0,
Subscript[c, "2223"] -> 0, Subscript[c, "2213"] -> 0, Subscript[c, "2212"] -> 0,
Subscript[c, "3323"] -> 0, Subscript[c, "3313"] -> 0, Subscript[c, "3312"] -> 0,
Subscript[c, "2311"] -> 0, Subscript[c, "2322"] -> 0, Subscript[c, "2333"] -> 0, Subscript[c, "2313"] -> 0, Subscript[c, "2312"] -> 0,
Subscript[c, "1311"] -> 0, Subscript[c, "1322"] -> 0, Subscript[c, "1333"] -> 0, Subscript[c, "1323"] -> 0, Subscript[c, "1312"] -> 0,
Subscript[c, "1211"] -> 0, Subscript[c, "1222"] -> 0, Subscript[c, "1233"] -> 0, Subscript[c, "1223"] -> 0, Subscript[c, "1213"] -> 0};

anisotropic // Normal // # /. orthotropicRules & // MatrixForm


This is very late to the party, but someone might still be interested. You can do what you want by solving the conditions of the material symmetry group. My answer is organized in

1) Background information

2) Mathematica code

1) Background information

First, some background information for those who are not engineers but want to understand better the question.

Index symmetries

A stiffness tensor $C$ is a fourth-order tensor with components $c_{ijkl}$ which maps symmetric second-order tensors into symmetric second-order tensors, i.e., $\sigma_{ij} = c_{ijkl} \varepsilon_{kl}$ (linear elastic law), $\sigma$ (stress) and $\varepsilon$ (strain) being arbitrary symmetric second-order tensors. Due to the symmetry of the second-order tensors, $C$ is allowed to be minor symmetric, i.e., $c_{ijkl} = c_{jikl} = c_{ijlk}$. The not minor symmetric part of $C$ is irrelevant for the elastic law and is dropped. If the stress $\sigma$ is to be related to an elastic energy potential $W$ (referred to as hyperelastic behavior), i.e., $\sigma = \partial W / \partial \varepsilon$, then, due to Schwarz's theorem, the stiffness tensor $c_{ijkl} = \partial^2 W / \partial \varepsilon_{ij} \partial \varepsilon_{kl}$ has to possess the major symmetry, i.e., $c_{ijkl} = c_{klij}$.

Material symmetry

A material with stiffness $C$ is said to possess the material symmetry group $G$ (e.g., triclinic, orthotropic, transversally isotropic, ...) if

$$C = Q \star C \qquad Q \in G$$

holds, where $Q$ are second-order tensors, referred to as symmetry transformations of $C$. The product $\hat{C} = Q \star C$ (referred to here as Rayleigh product) is defined in components as

$$\hat{c}_{ijkl} = Q_{im}Q_{jn}Q_{ko}Q_{lp}c_{mnop}$$

For solids, $G$ is a subset of the orthogonal group. In solid mechanics, if suffices to consider rotation matrices $Q$ from the rotational group $SO(3)$. If $G = \{I\}$, $I$ being the identity matrix, then $C$ is said to triclinic. If $G$ possesses more than the identity transformation, then different material classes can be defined (different anisotropy types). If $G = SO(3)$, the $C$ is said to be isotropic (no direction dependency).

2) Mathematica code

Index symmetries

You might want to derive the symmetrized tensor from scratch. In the following code I will create a fourth-order tensor C0 and ask for its free components if I impose minor and major symmetries. The solution will be inserted into C0 and I will check the number of linearly independent variables / degrees of freedom dof, which is 21 in this case (as also known in the mechanics community).

C0 = Normal@
SymmetrizedArray[
pos_ :> Subscript[c, StringJoin[ToString /@ pos]], {3, 3, 3,
3}, {}];
csol = Solve[
{
C0 == TensorTranspose[C0, {2, 1, 3, 4}](*minor symmetries*)
, C0 == TensorTranspose[C0, {1, 2, 4, 3}]
, C0 == TensorTranspose[C0, {3, 4, 1, 2}](*major symmetry*)
}
][[1]];
C0 = C0 /. csol;
dof = Variables@C0;
Length@dof


21

If you are familiar with the normalized Voigt notation, see, e.g., this paper, you can represent a minor symmetric fourth-order tensor as a 6x6 matrix (I will only describe this more in detail if needed, since this is just an extra not necessarily asked by the OP)

$$\left( \begin{array}{cccccc} c_{1111} & c_{1122} & c_{1133} & \sqrt{2} c_{1123} & \sqrt{2} c_{1113} & \sqrt{2} c_{1112} \\ c_{2211} & c_{2222} & c_{2233} & \sqrt{2} c_{2223} & \sqrt{2} c_{2213} & \sqrt{2} c_{2212} \\ c_{3311} & c_{3322} & c_{3333} & \sqrt{2} c_{3323} & \sqrt{2} c_{3313} & \sqrt{2} c_{3312} \\ \sqrt{2} c_{2311} & \sqrt{2} c_{2322} & \sqrt{2} c_{2333} & 2 c_{2323} & 2 c_{2313} & 2 c_{2312} \\ \sqrt{2} c_{1311} & \sqrt{2} c_{1322} & \sqrt{2} c_{1333} & 2 c_{1323} & 2 c_{1313} & 2 c_{1312} \\ \sqrt{2} c_{1211} & \sqrt{2} c_{1222} & \sqrt{2} c_{1233} & 2 c_{1223} & 2 c_{1213} & 2 c_{1212} \\ \end{array} \right)$$

Please note, that the matrix just given is only for a minor symmetric fourth-order tensor, for additionally major symmetry the 6x6 matrix is also symmetric (to be given in a moment). I still wanted to give this matrix representation, since it will be useful for the material symmetries and is used in almost all standard mechanics text books.

The following code converts a fourth-order tensor into the normalized Voigt notation

nvn4[A_] := Module[
{ip, i, j},
ip = {{1, 1}, {2, 2}, {3, 3}, {2, 3}, {1, 3}, {1, 2}};
Table[
A[[
ip[[i, 1]]
, ip[[i, 2]]
, ip[[j, 1]]
, ip[[j, 2]]
]]*If[4 <= i <= 6, Sqrt[2], 1]*If[4 <= j <= 6, Sqrt[2], 1]
, {i, 6}, {j, 6}]
];
nvn4::usage =
"nvn4[A] returns the normalized Voigt notation of the fourth-order \
tensor A.";


You can then use this on C0

nvn4[C0] // MatrixForm


where you can see, that the 6x6 matrix is symmetric (due to the imposed major symmetry of the fourth-order tensor).

Material symmetry

I will need the Rayleigh product

rp[A_, B_] := Block[
{n, it, t1},
n = TensorRank[B];
it = RotateLeft@Range[n];
t1 = B;
Do[t1 = TensorTranspose[A.t1, it], {i, n}];
t1
];
rp::usage =
"rp[A,B] computes the Rayleigh product of the tensor B with \
transformation tensor A.";


Let's define the symmetry group for orthotropic tensors

G = {
RotationMatrix[Pi, {1, 0, 0}]
, RotationMatrix[Pi, {0, 1, 0}]
, RotationMatrix[Pi, {0, 0, 1}]
};


You can then just use Solve based on the material symmetry conditions, insert the result into a local stiffness Cloc, take a look at the number of degrees of freedom and at the 6x6 matrix.

Gsol = Solve[Table[rp[G[[i]], C0] == C0, {i, Length@G}]][[1]];
Cloc = C0 /. Gsol;
dof = Variables@Cloc;
Length@dof
Cloc // nvn4 // Simplify // MatrixForm


Now you know that linear elastic orthotropic materials have 9 free components and can be represented by the just given matrix.

Naturally, you just have to know the symmetry transformations for the case you are taking a look at, e.g., for cubic materials you will need to know the 24 symmetry transformations.

• in some books, they write C_2211, C_3311, C_1211 C_2311, ... Jul 12, 2018 at 10:04
• Yes, but since you have the major symmetry c_ijkl = c_klij, you get the same results, e.g., c_3311 = c_1133 and so on. Jul 12, 2018 at 10:27
• okay I found that C_2211, C_3311, C_1211 C_2311 are used in many many books... this paper I cannot open, can you give us a new links? thanking you! Jul 12, 2018 at 10:29
• Right now I only have my phone with me. I will give you another link to a paper today evening (Germany). Jul 12, 2018 at 10:52
• or you can tell me the name of this paper, I can download via VPN. Jul 12, 2018 at 11:06