The article under concern is:
Germano, M. (1992). Turbulence: The filtering approach. Journal of Fluid Mechanics, 238, 325-336. doi:10.1017/S0022112092001733
On page 328, Germano considers the Navier-Stokes equation (equation 17) written in Einstein notation:
$u_{i, t}+\left(u_{i} u_{k}\right)_{, k}=-p_{, i}+\sigma_{i k, k}$
and states that 'by taking a moment of (17) with $u_j$ and adding this to another moment of the same equation but with the indices interchanged', the following equation can be achieved:
$\left(u_{i} u_{j}\right)_{, t}+\left(u_{i} u_{j} u_{k}\right)_{, k}=-\left[p u_{i} \delta_{j k}+p u_{j} \delta_{i k}-\nu\left(u_{i} u_{j}\right)_{, k}\right]_{,k}+2 p s_{i j}-2 v u_{i, k} u_{j, k}$
where $\sigma_{i j}=2 \nu s_{i j} ; \quad s_{i j}=\frac{1}{2}\left(u_{i, j}+u_{j, i}\right)$.
The derivation is made with the following commuting properties of the moment operator: $\left( f_{, t}\right)=(f)_{, t} ; \quad\left( f_{, k}\right)=( f)_{, k}$
After hours of attempts, I got nowhere near the final derivation. One of the significant issues I am facing is how $p$, without the derivative (so $p_{,i}$), is achieved as well as where the $\delta_{jk}$ and $\delta_{ik}$ come from. Also, regarding the 'indices interchanged' part, is the author referring to equation (17) itself or after taking the moment of it with $u_j$?
I would really appreciate if someone could shed some light on how I should go about deriving this.
I'm going to do this not using Einstein summation and with slightly different notation for partial derivatives. I have a feeling the notation is part of the difficulty. In my notation the NSE will take the form
$$\partial_t u_i + \sum_{k}\partial_{k}(u_i u_k)=-\partial_i p+\sum_{k}\partial_k\sigma_{ik}\ .$$
Now multiply by $u_j$ and write the same formula with $j\leftrightarrow i$. This is
$$\begin{align}u_j\partial_t u_i + \sum_{k}u_j\partial_{k}(u_i u_k)&=-u_j\partial_i p+\sum_{k} u_j\partial_k\sigma_{ik}\ ,\\ u_i\partial_t u_j + \sum_{k}u_i\partial_{k}(u_j u_k)&=-u_i\partial_j p+\sum_{k} u_i\partial_k\sigma_{jk}\ . \end{align}$$
Now adding these and using the fact that $\partial_t(fg)=g\partial_tf+f\partial_tg$ we obtain
$$\partial_t(u_i u_j) +\sum_{k} \left[u_j\partial_k(u_iu_k)+u_i\partial_k(u_ju_k)\right]=-u_j\partial_i p-u_i\partial_jp+\sum_{k}\left[u_j\partial_k\sigma_{ik}+u_i\partial_k\sigma_{jk} \right]\ .$$
The $2$nd term on the left hand side is computed using the divergence free condition $\sum_{k}\partial_ku_k=0$. With this in mind
$$\begin{align}\sum_{k}\partial_k(u_i u_j u_k) &= \sum_{k} u_i\partial_k(u_j u_k)+\sum_{k} u_j u_k\partial_k u_i\\ &= \sum_{k} u_i\partial_k(u_j u_k)+\sum_{k} u_j u_k\partial_k u_i+\sum_{k}u_i u_j\partial_k u_k\\ &=\sum_{k} u_i\partial_k(u_j u_k)+\sum_{k} u_j\partial_k(u_iu_k)\ . \end{align}$$
On the right hand side
$$\begin{align}-u_j\partial_ip-u_i\partial_j p &= -\partial_i(pu_j)+p\partial_iu_j-\partial_j(p u_i)+p\partial_ju_i\\ &=-\partial_i(pu_j)-\partial_j(p u_i)+2ps_{ij}\\ &=-\left[\sum_{k}\partial_k(pu_j\delta_{ik})+\sum_{k}\partial_k(pu_i\delta_{jk})\right]+2p s_{ij}\ .\end{align}$$
Finally we need to compute the last term on the right hand side
$$\begin{align}\sum_{k}\left[u_j\partial_k\sigma_{ik}+u_i\partial_k\sigma_{jk}\right]&=\nu\sum_{k}\left[u_j\partial_k(\partial_ku_i+\partial_i u_k)+u_i\partial_k(\partial_k u_j+\partial_j u_k)\right]\\ &=\nu\sum_{k}\left[u_j\partial_k^2 u_i+u_i\partial_k^2 u_j\right]\ ,\end{align}$$
by Clairaut's theorem and the divergence free condition. We continue
$$\begin{align}\nu\sum_{k}\left[u_j\partial_k^2 u_i+u_i\partial_k^2 u_j\right]=\nu\left[\sum_{k}\partial_k^2(u_i u_j)-2\sum_{k}\partial_k u_i\partial_k u_j\right]\ .\end{align}$$
Putting this together gives
$$\begin{align}\partial_t(u_iu_j)+\sum_{k}\partial_k(u_i u_j u_k)=&-\left[\sum_{k}\{\partial_k(pu_i\delta_{jk})+\partial_k(pu_j\delta_{ik})-\nu\partial_k^2(u_iu_j)\right]\\&+2ps_{ij}-2\nu\sum_{k}\partial_k u_i \partial_k u_j\end{align}$$