Solution (source code)

= Solution

Using the given linear momentum balance, multiply by the <velocity> and integrate over the fixed volume:
$$
\frac{d}{dt}\int_V\frac12\rho u_iu_i\,dV=\int_Vu_i\partial_j\sigma_{ij}\,dV
=\int_{\partial V}u_i\sigma_{ij}n_j\,dS-\int_V\sigma_{ij}\partial_j u_i\,dV.
$$
This is <integration by parts> component by component, or the <divergence theorem>. The first term is boundary <stress> power. With component-first <velocity gradient> $K_{ij}=\partial_j u_i$, the second is the <double contraction> $\Phi=\sigma:K$.

For a <Newtonian fluid>, $\sigma=-pI+2\mu e$, where $e=(K+K^T)/2$ is the <rate-of-strain tensor> and $\mu\ge0$ is <dynamic viscosity>. <Incompressibility> removes the pressure contribution, and the <symmetric tensor> $e$ is orthogonal to the antisymmetric part of $K$. Therefore
$$
\boxed{\Phi=2\mu e:e=2\mu\sum_{ij}e_{ij}^2\ge0.}
$$
The displayed energy identity follows for the stated linear balance. A full advective momentum equation in a fixed control volume would also contribute an advective <kinetic energy> flux. The nonnegativity statement here is for a purely Newtonian viscous stress; a <viscoelastic> stress can exchange stored <free energy> with the flow.