You should start with the strain energy density $\psi$, then define:
$$
C_{ijkl} = \frac{\partial^2 \psi}{\partial \epsilon_{ij}\partial \epsilon_{kl}}, $$
and then define $$
\sigma_{ij} = C_{ijkl} \epsilon_{kl}
$$
The remainder of my answer will be about explaining why you have to do it that way. Firstly it is physical, there really is energy associated with strain, and if there weren't there would not be any stress. Secondly, it is linear exactly because we are considering the changes in energy due to small strains.
But let's go back to the minor symmetries. We need $C_{ijkl}=C_{jikl}$ because otherwise $\sigma_{ij}\neq\sigma_{ji}$ (and then we get arbitrarily large, hence unphysical, angular velocities for smaller and smaller regions). But the other minor symmetry is not required. If someone handed you a random rank four tensor, let's call it $B_{ijkl}$, and called it an elastic tensor and it didn't have the second minor symmetry, you can define $C_{ijkl}=(B_{ijkl}+B_{ijlk})/2$ $D_{ijkl}=(B_{ijkl}-B_{ijlk})/2$ and then $B=C+D$ but when D is contracted with a symmetric rank two tensor (like the strain tensor) it gives zero. So the part of the rank four tensor without that second minor symmetry simply doesn't contribute, as an actor it does nothing (when you think all the elasticity does is give you stress from strain). So you might as well assume your tensor has both minor symmetries because it acts like it has the second ($B$ and $C$ act the same on symmetric tensors) and it has to have the first.
Did I bring that up to be pedantic? No, I brought it up because the same thing happens if you contract the elasticity tensor with a rank four symmetric combination of the strain tensor. The part of the elasticity tensor without the major symmetry doesn't contribute to the strain energy density. So a random tensor needs the first minor symmetry to be physical. But you might as well assume it has the second minor symmetry since it doesn't affect the stress-strain relationship. And you might as well assume it has the major symmetry because the part that doesn't will not contribute to the strain energy density.
But it is the strain energy density that is physical, and how it changes is what elasticity is. So you aren't really deriving these symmetries as much as saying that only the symmetric ones generate the physical things you want, energy when given strain. And a real derivation should start with strain energy density and strain, and then just define elasticity from that.
First, I'm going to define some notation. I prefer the direct notation of modern continuum mechanics as presented in this book by Gurtin et al., so I won't be carrying around the indices that you use. Lastly, the following holds for small deformation, rate-independent plasticity with isotropic hardening.
- $\mathbf{T}$ = Cauchy stress tensor
- $\mathbf{T}_0 = \mathbf{T} - \frac{1}{3}\mathrm{tr}(\mathbf{T}) \mathbf{1}$ = Deviatoric part of Cauchy stress
- $\bar{\sigma} = \sqrt{\frac{2}{3} \mathbf{T}_0 : \mathbf{T}_0}$ = Equivalent tensile stress
- $\mathbf{N}^p = \sqrt{\frac{3}{2}}\frac{\mathbf{T}_0}{\bar{\sigma}}$ = Direction of plastic flow (assume co-directionality)
- $\mathbf{E} = \mathbf{E}^e + \mathbf{E}^p$ = Additive strain decomposition
- $\mathbf{T} = \mathbb{C}:\mathbf{E}^e = (2G)\mathbf{E}^e + \lambda\mathrm{tr}(\mathbf{E}^e)\mathbf{1}$ = Linearized Elasticity constitutive equation
- $\dot{\mathbf{E}}^p = \sqrt{\frac{3}{2}}\dot{\bar{\epsilon}}^p \mathbf{N}^p$
- $\dot{\bar{\epsilon}}^p = \sqrt{\frac{2}{3}} |\dot{\mathbf{E}}^p|$ = equivalent plastic strain rate
- $\bar{\epsilon}^p = \int_0^t \dot{\bar{\epsilon}}^p(t')\,dt'$ = accumulated plastic strain
- $Y(\bar{\epsilon}^p)$ = "flow resistance" = 1D tensile yield stress $\sigma_y$ (function of $\bar{\epsilon}^p$)
- $H(\bar{\epsilon}^p) = \frac{dY}{d\bar{\epsilon}^p}$ = strain-hardening rate
- $f(\mathbf{T}_0, Y(\bar{\epsilon}^p)) = \bar{\sigma} - Y(\bar{\epsilon}^p) \le 0 \; \forall \,\mathbf{T}, \,Y$ <-- Yield Function
Next, we write out the Kuhn-Tucker consistency condition at yield (which occurs at $f=0$):
- If, $f = 0$ then $\dot{\bar{\epsilon}}^p \dot{f} = 0$.
This may seem like a strong statement, and it is, but if you think about all of the possible types of loading from a point on the "yield surface" (defined as all stress states that give $f=0$), you will find this to be true.
Now, we apply this condition to find $\dot{\bar{\epsilon}}^p$.
First, compute $\dot{f}$.
\begin{align}
\dot{f} &= \sqrt{\frac{3}{2}} \dot{\bar{\sigma}} - \dot{Y}\\
&= \sqrt{\frac{3}{2}} \frac{\dot{\mathbf{T}}_0 : \mathbf{T}_0}{|\mathbf{T}_0|} - H(\bar{\epsilon}^p)\dot{\bar{\epsilon}}^p\\
&= \sqrt{\frac{3}{2}} \dot{\mathbf{T}}_0 : \mathbf{N}^p - H(\bar{\epsilon}^p)\dot{\bar{\epsilon}}^p \\
&= \left(\sqrt{\frac{3}{2}}\right)(2G)[\dot{\mathbf{E}} - \dot{\mathbf{E}}^p] : \mathbf{N}^p - H(\bar{\epsilon}^p)\dot{\bar{\epsilon}}^p \\
&= \left(\sqrt{\frac{3}{2}}\right)(2G)\dot{\mathbf{E}} : \mathbf{N}^p - (3G + H)\dot{\bar{\epsilon}}^p
\end{align}
Now, we have to consider three cases: 1) Elastic unloading, 2) Neutral Loading, 3) Plastic Loading
- Elastic Unloading: This is like releasing the stress on a body that is at its yield point. Thus, $\dot{\bar{\epsilon}}^p = 0$, and $\dot{\mathbf{E}}:\mathbf{N}^p < 0$. Therefore, $\dot{f} < 0$.
- Neutral Loading: This one is a bit trickier. Here, the loading is tangent to the yield surface. For this case, the strain rate is orthogonal to the direction of plastic flow, thus $\dot{\mathbf{E}}:\mathbf{N}^p = 0$. Here, still, we have not loaded the body plastically, so $\dot{\bar{\epsilon}}^p = 0$. In codes, this is usually handled without special consideration by either the above or below cases, with the same result.
- Plastic Loading: Here, we are loading the body plastically. This means that $\dot{\mathbf{E}}:\mathbf{N}^p > 0$. Since $f=0$ at the yield surface, and $f$ cannot have a positive value, this means that $\dot{f} = 0$. This is how we can find the value of $\dot{\bar{\epsilon}}^p$.
$$
\dot{\bar{\epsilon}}^p = \frac{\sqrt{3/2} (2G) (\dot{\mathbf{E}}:\mathbf{N}^p)}{3G + H(\bar{\epsilon}^p)}
$$
I glossed over some of the more tedius algebra and buried the assumption that $(3G + H)>0$, which is true for every case that I have encountered. You can go on from here to compute the so-called "elasto-plastic" modulus, which is necessary to do implicit integration in a finite element method.
Best Answer
I use this notation
the transformation matrix, transformed a vector components from rotate system index B to inertial system index I
rotation about the x-axis angle $~\alpha~$ between y and y'
$${_B^I}\,Q_x=\left[ \begin {array}{ccc} 1&0&0\\ 0&\cos \left( \alpha \right) &-\sin \left( \alpha \right) \\ 0& \sin \left( \alpha \right) &\cos \left( \alpha \right) \end {array} \right] $$
rotation about the y-axis angle $~\alpha~$ between x and x'
$${_B^I}\,Q_y= \left[ \begin {array}{ccc} \cos \left( \alpha \right) &0&\sin \left( \alpha \right) \\ 0&1&0\\ -\sin \left( \alpha \right) &0&\cos \left( \alpha \right) \end {array} \right] $$
rotation about the z-axis angle $~\alpha~$ between x and x'
$${_B^I}\,Q_z=\left[ \begin {array}{ccc} \cos \left( \alpha \right) &-\sin \left( \alpha \right) &0\\ \sin \left( \alpha \right) &\cos \left( \alpha \right) &0\\ 0&0&1\end {array} \right] $$
vector transformation from B to I system
$$\mathbf v_I={_B^I}\mathbf Q\,\mathbf v_B$$
matrix transformation
$$\mathbf M_I=\mathbf ={_B^I}\mathbf Q \,\mathbf M_B\,{_I^B}\mathbf Q =\mathbf Q\,\mathbf M_B\,\mathbf Q^T\\ \mathbf M_B=\mathbf ={_I^B}\mathbf Q \,\mathbf M_I\,{_B^I}\mathbf Q =\mathbf Q^T\,\mathbf M_I\,\mathbf Q $$
your matrix is
$$\mathbf \epsilon_I= \left[ \begin {array}{ccc} \epsilon_{{11}}&0&0\\ 0& \epsilon_{{22}}&0\\ 0&0&\epsilon_{{22}}\end {array} \right] \\ \mathbf \epsilon_B=Q^T\,\mathbf \epsilon_I\,\mathbf Q $$
for $~\mathbf Q=\mathbf Q_x~$ you obtain
$$\mathbf \epsilon_B=\mathbf \epsilon_I$$
for $~\mathbf Q=\mathbf Q_y~$
$$\mathbf \epsilon_B=\left[ \begin {array}{ccc} \left( \cos \left( \alpha \right) \right) ^{2}\epsilon_{{11}}+\epsilon_{{22}}- \left( \cos \left( \alpha \right) \right) ^{2}\epsilon_{{22}}&0&\cos \left( \alpha \right) \sin \left( \alpha \right) \left( -\epsilon_{{22}}+\epsilon_ {{11}} \right) \\ 0&\epsilon_{{22}}&0 \\ \cos \left( \alpha \right) \sin \left( \alpha \right) \left( -\epsilon_{{22}}+\epsilon_{{11}} \right) &0& \left( \cos \left( \alpha \right) \right) ^{2}\epsilon_{{22}}+\epsilon_{{11} }- \left( \cos \left( \alpha \right) \right) ^{2}\epsilon_{{11}} \end {array} \right] $$
for $~\mathbf Q=\mathbf Q_z~$
$$\mathbf \epsilon_B= \left[ \begin {array}{ccc} \left( \cos \left( \alpha \right) \right) ^{2}\epsilon_{{11}}+\epsilon_{{22}}- \left( \cos \left( \alpha \right) \right) ^{2}\epsilon_{{22}}&-\cos \left( \alpha \right) \sin \left( \alpha \right) \left( -\epsilon_{{22}}+\epsilon_ {{11}} \right) &0\\ -\cos \left( \alpha \right) \sin \left( \alpha \right) \left( -\epsilon_{{22}}+\epsilon_{{11}} \right) & \left( \cos \left( \alpha \right) \right) ^{2}\epsilon_{{ 22}}+\epsilon_{{11}}- \left( \cos \left( \alpha \right) \right) ^{2} \epsilon_{{11}}&0\\ 0&0&\epsilon_{{22}}\end {array} \right] $$