Assume that K is a diffusion coefficient for momentum diffusion.
Assume we are in two-dimensional space (horizontal divergence and vertical vorticity are relevant) on the horizontal plane. Then the friction term in the momentum equation reads either
$$\nabla\cdot(K(\nabla\mathbf{v}+\mathbf{v}\nabla)-K\nabla\cdot\mathbf{v}\mathbb{E})$$
or, after some algebra
$$\nabla(K\nabla\cdot\mathbf{v})-\nabla\times(K(\nabla\times\mathbf{v})-2(\nabla K)\nabla\cdot\mathbf{v}+2\nabla K\cdot \nabla \mathbf{v}+2\nabla K\times(\nabla\times\mathbf{v})$$
In more detail:
\begin{eqnarray}
\mathbf{F}&=&\nabla\cdot(K(\nabla\mathbf{v}+\mathbf{v}\nabla)
-K\nabla\cdot\mathbf{v}\mathbb{E})\nonumber\\
&=&\nabla K\cdot\nabla\mathbf{v}+K\nabla\cdot\nabla\mathbf{v}+\nabla K\cdot\mathbf{v}\nabla+K\nabla\cdot\mathbf{v}\nabla-\nabla(K\nabla\cdot\mathbf{v})\nonumber\\
&=&\nabla\mathbf{v}\cdot\nabla K-\nabla K\times(\nabla\times\mathbf{v})
+K\nabla(\nabla\cdot\mathbf{v})-K\nabla\times(\nabla\times\mathbf{v})
\nonumber\\&&+\nabla\mathbf{v}\cdot\nabla K +K \nabla\nabla\cdot\mathbf{v}-\nabla(K\nabla\cdot\mathbf{v})\nonumber\\
&=&-\nabla\times (K(\nabla\times\mathbf{v}))+2\nabla\mathbf{v}\cdot\nabla K
+2K\nabla\nabla\cdot\mathbf{v}-\nabla(K\nabla\cdot\mathbf{v})\nonumber\\
&=&-\nabla\times (K(\nabla\times\mathbf{v}))+2\nabla\mathbf{v}\cdot\nabla K
-\nabla(K\nabla\cdot\mathbf{v})+\\&&+2\nabla(K\nabla\cdot\mathbf{v})-2(\nabla K)\nabla\cdot\mathbf{v}\nonumber\\
&=&-\nabla\times (K(\nabla\times\mathbf{v}))+\nabla(K\nabla\cdot\mathbf{v})
+2\nabla\mathbf{v}\cdot\nabla K-2(\nabla K)\nabla\cdot\mathbf{v}\nonumber\\
&=&-\nabla\times (K(\nabla\times\mathbf{v}))+\nabla(K\nabla\cdot\mathbf{v})
-2(\nabla K)\nabla\cdot\mathbf{v}+\\&&+2\nabla K\cdot\nabla\mathbf{v}+2\nabla K\times(\nabla\times\mathbf{v})
\end{eqnarray}
This has been checked by another person meanwhile.
I have two problems with those equivalent forms.
First, when computing the energy dissipation for a constant diffusion coefficient, I get two different answers, depending on the chosen form, namely
either $K(D^2+\zeta^2)$ where $D$ is the divergence and $\zeta$ is the vorticity
or $K(E^2+F^2)$ where $F=\partial_xv+\partial_yu$ and $E=\partial_xu-\partial_yv$ are the shear and strain deformations.
Which is the right dissipation and why?
Second, the form with $D$ and $\zeta$ is 'vector invariant' or, better said, it only relies on Gauss and Stokes theorem, and is independent of the chosen coordinate system. This seems not to be the case for the form with the deformations: I personally use in my modeling not a cartesian, but a trivariate coordinate system (my grid boxes are hexagons). Then, the vector components are linearly dependent. I stumbled across the problem, that the numerical discretisation gives 2 times the expected Laplacian when I use the shear and strain deformations in the tensor, but it gives the right Laplacian if I use the divergence and vorticity in the tensor.
What is going on here?