Tensor deviator calculation rules In the field of continuum/solid mechanics, there are often deviatoric tensors defined, like for the derivation (comma in einstein notation) of a displacement
$$\mathrm{dev}(u_{i,j})=u_{i,j}-\frac{1}{3}\,u_{k,k}\,\delta_{ij},$$ where $\delta_{ij}$ is the Kronecker delta.
If one performs a double contraction with another tensor, for example $\sigma_{ij}$, I observed, that several textbooks move the deviator operator to the other tensor, namely
$$\sigma_{ij}\,\mathrm{dev}(u_{i,j}) = \mathrm{dev}(\sigma_{ij})\,u_{i,j} = \mathrm{dev}(\sigma_{ij})\,\mathrm{dev}(u_{i,j})$$
I (engineering student ;) ) checked this relation with Python (see code below), but I did not find any mathematical explaination therefore in literature. Can you give me a hint in which books I have to look in oder to find calculation rules like for the manipulation of terms with deviatoric tensors?
The abovementioned Python sourcecode:
# -*- coding: utf-8 -*-
import numpy as np   

def dev(tens):
    return tens - 1./len(tens)*np.trace(tens)*np.eye(len(tens))

u = np.random.rand(3,3)
sigma = np.random.rand(3,3)

print("u * dev sigma")
print(np.einsum('ij,ij',dev(sigma),u))
print("\ndev u * sigma")
print(np.einsum('ij,ij',sigma,dev(u)))
print("\ndev u * dev sigma")
print(np.einsum('ij,ij',dev(sigma),dev(u)))
print("\n u * sigma")
print(np.einsum('ij,ij',sigma,u))


$\require{cancel}$
First note that
$$
\delta_{ii} = n \tag{1}
$$
where $n = 3$ is the dimension of your tensors. With this in mind
\begin{eqnarray}
\sigma_{ij}~{\rm dev}(u_{ij}) &=& \sigma_{ij} \left(u_{ij} - \frac{1}{n}u_{kk}\delta_{ij} \right) \\
&=& \sigma_{ij}u_{ij} - \frac{1}{n} \sigma_{ij}u_{kk} \delta_{ij} \\
&=& \sigma_{ij}u_{ij} - \frac{1}{n} \sigma_{ii} u_{kk}  \tag{2}
\end{eqnarray}
And
\begin{eqnarray}
{\rm dev}(\sigma_{ij})~{\rm dev}(u_{ij}) &=& \left(\sigma_{ij} - \frac{1}{n}\sigma_{ll}\delta_{ij} \right) \left(u_{ij} - \frac{1}{n}u_{kk}\delta_{ij} \right) \\
&=& \sigma_{ij}u_{ij} - \frac{1}{n}\sigma_{ii} u_{kk} - \frac{1}{n}\sigma_{ll} u_{ii} + \frac{1}{n^2}\sigma_{ll} u_{kk} \delta_{ii} \\
&\stackrel{(1)}{=}& \sigma_{ij}u_{ij} - \frac{1}{n}\sigma_{ii} u_{kk} - \cancel{\frac{1}{n}\sigma_{ll} u_{ii}} + \cancel{\frac{1}{n^2}\sigma_{ll} u_{kk} n} \\
&\stackrel{(2)}{=}& \sigma_{ij}~{\rm dev}(u_{ij}) \tag{3}
\end{eqnarray}
