As you suspect, the mechanical energy equation you derive is missing terms that account for viscous effects. The other terms you derive are correct.
The RHS should be $\nabla \cdot \mathbf{\sigma}$, the divergence of the stress tensor, and you have edited accordingly.
This clearly is the correct term to include because the rate of increase of momentum in a fixed volume $V$ of fluid should be related to the viscous stress force of the surrounding fluid acting on the bounding surface of the volume $\partial V$,
$$\text{rate of momentum change } = \int_{\partial V}\mathbf{\sigma} \cdot \mathbf{n} \, dS = \int_V \nabla \cdot \mathbf{\sigma} \, dV.$$
Your error stems from incorrectly forming the dot product of the velocity and the viscous term $\eta\,\nabla \cdot [\nabla \mathbf{v} + \nabla \mathbf{v}^T] = \eta \, \nabla^2 \mathbf{v} + \eta \, \nabla (\nabla \cdot \mathbf{v})$.
In particular, you are missing $ \mathbf{v} \cdot [\eta\,\nabla \cdot \nabla \mathbf{v}]$ which should give rise to the viscous dissipation term, representing the action of viscous forces in converting useful mechanical energy into heat.
The correct treatment of that term is
$$\tag{*} \eta\,\mathbf{v} \cdot\,\nabla \cdot [\nabla \mathbf{v} + \nabla \mathbf{v}^T] = \underbrace{-\eta \nabla \mathbf{v} \,\mathbf{:}\, \nabla \mathbf{v}}_{\text{viscous dissipation}} + \frac{\eta}{2}\nabla^2(|\mathbf{v}|^2) + \eta \, \mathbf{v} \cdot \nabla (\nabla \cdot \mathbf{v}), $$
where the double dot product represents $\sum_{i=1}^3 \sum_{j=1}^3 |\frac{\partial v_i}{\partial x_j}|^2$, a norm of the velocity gradient, and $|\mathbf{v}|$ is the usual Euclidean norm of the velocity vector.
For an easy derivation, we can expand $\mathbf{v} \cdot\,\nabla \cdot [\nabla \mathbf{v}+ \nabla \mathbf{v}^T] $ using Cartesian coordinates and the Einstein convention (repeated index implies summation over that index). First we obtain
$$\nabla \cdot [\nabla \mathbf{v}+ \nabla \mathbf{v}^T]_{i} = \partial_j(\partial_j v_i + \partial_iv _j) = \partial_j(\partial_j v_i) + \partial_i(\partial_jv_j),$$
and forming the dot product we get
$$ \mathbf{v} \cdot\,\nabla \cdot [\nabla \mathbf{v} + \nabla \mathbf{v}^T] = v_i\,\partial_j(\partial_j v_i) + v_i \,\partial_i(\partial_jv_j) \\ = -\partial_jv_i \, \partial_j v_i + \partial_j[v_i(\partial_j v_i)] + v_i \,\partial_i(\partial_jv_j) \\ = -\partial_jv_i \, \partial_j v_i + \frac{1}{2}\partial_j\,\partial_j (v_iv_i) + v_i \,\partial_i(\partial_jv_j), $$
which is the component form of the coordinate-free term on the RHS of (*).
The most generic form of the Navier-Stokes equation is
$$
\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot(\rho \mathbf{u} \otimes \mathbf{u}) = \nabla p - \mathbf{f} + \nabla \cdot \mathbf{S},
\tag{1}
$$
in which $\mathbf{S}$ is the shear stress tensor ($\mathbf{S}=\mu\nabla \mathbf{u}$ in your case). The continuity equation is
$$
\frac{\partial \rho}{\partial t}+\nabla \cdot (\rho \mathbf{u})=0.
\tag{2}
$$
Using the product rule for derivatives in equation $(1)$ leads to
$$
\rho \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u}\frac{\partial \rho}{\partial t} + \rho \mathbf{u} \cdot \nabla \mathbf{u} + \mathbf{u} \nabla \cdot(\rho \mathbf{u}) = \nabla p - \mathbf{f} + \nabla \cdot \mathbf{S},
\tag{3}
$$
and using the continuity equation we see that the second and fourth terms in LHS cancel each other, leading to
$$
\rho\left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = \nabla p - \mathbf{f} + \nabla \cdot \mathbf{S}.
\tag{4}
$$
Since it's usually assumed that the continuity equation holds, equation $(4)$ is completly equivalent to equation $(1)$. In practice, however, there are some "differences":
Equation $(1)$ is called the conservative form of Navier-Stokes equation, while equation $(4)$ is called the non-conservative form. These names are a bit misleading: both equations are the momentum conservation equations. However, when solving numerically the governing equations of fluid dynamics, it's sometimes more useful to use equation $(1)$. It's basically due to the fact that accross a shock wave the velocity $\mathbf{u}$ is discontinuous (and, therefore, equation $(4)$ has the gradient of a discontinuous function), while $\rho \mathbf{u} \otimes \mathbf{u}$ is continuous even accross the shock. See that equation $(1)$ is called conservative form because it has derivatives of a conserved quantity (the momentum, i.e., $\rho \mathbf{u}$) while equation $(4)$ has derivatives of a non-conserved quantity (the velocity).
Equation $(4)$ explicitly shows the transport of momentum in the term $\mathbf{u} \cdot \nabla \mathbf{u}$. Notice that the conservation equation of the property $\phi$ (which can be enthalpy, vorticity, chemical species, etc.) will have a term $\mathbf{u} \cdot \nabla \phi$. Therefore, equation $(4)$ "looks like" every other conservation equation. It's usually defined the material derivative of property $\phi$ as
$$
\frac{D \phi}{Dt} = \frac{\partial\phi}{\partial t} + \mathbf{u} \cdot \nabla\phi,
$$
which can be interpreted as the rate of change of the property $\phi$ along time in a particle of fluid while this particle is transported by the flow. Then, the conservation equation of any property $\phi$ can be written generically as
$$
\frac{D \phi}{Dt} = \text{source terms},
$$
in which the source terms for the case $\phi=\mathbf{u}$ (i.e., the Navier-Stokes equation) are $(\nabla p - \mathbf{f} + \nabla \cdot \mathbf{S})/\rho$.
Summarizing: both equations are the same. When you need to solve them numerically, equation $(1)$ can be more suitable. If you want to interpret the physical meaning of the terms of Navier-Stokes equation, equation $(4)$ is more suitable.
Best Answer
First, you need to recognize that $W_t$ represents a material control volume that is time-dependent as it moves and deforms with the fluid. Then you pass the derivative under the integral on the right-hand side by applying Reynolds transport theorem to obtain,
$$\frac{d}{dt} \int_{W_t} \rho \mathbf{u} \, dV = \int_{W_t} \frac{\partial (\rho \mathbf{u})}{\partial t} \, dV + \int_{\partial W_t}\rho \mathbf{u} \mathbf{u} \cdot \mathbf{n} \, dS,$$
where $\mathbf{n}$ is the outer unit normal vector field at the surface $\partial W_t$ bounding the control volume and the integral on the far right-hand side is a surface integral.
By the divergence theorem, it follows that
$$\tag{1}\frac{d}{dt} \int_{W_t} \rho \mathbf{u} \, dV = \int_{W_t} \frac{\partial (\rho \mathbf{u})}{\partial t} \, dV + \int_{W_t}\nabla \cdot(\rho \mathbf{u} \mathbf{u}) \,dV = \int_{W_t} \left[\frac{\partial (\rho \mathbf{u})}{\partial t} +\nabla \cdot(\rho \mathbf{u} \mathbf{u})\right] \,dV $$
The right-hand side you have written is incorrect. This should represent the net force acting on the control volume due to pressure and stress at the surface $\partial W_t$. Correctly written it is
$$\int_{\partial W_t}(-p\mathbf{n} + \mathbf{\sigma} \cdot \mathbf{n}) \, dS,$$
where $p$ is the pressure field and $\sigma $ is the (deviatoric) stress tensor. Applying the divergence theorem, we get
$$\tag{2}\int_{\partial W_t}(-p\mathbf{n} + \mathbf{\sigma} \cdot \mathbf{n}) \, dS = \int_{W_t} (- \nabla p + \nabla \cdot \sigma) \, dV$$
Equating (1) and (2), we have
$$\int_{W_t} \left[\frac{\partial (\rho \mathbf{u})}{\partial t} +\nabla \cdot(\rho \mathbf{u} \mathbf{u})\right] \,dV = \int_{W_t} (- \nabla p + \nabla \cdot \sigma) \, dV$$
Since this is true for any control volume $W_t$ we can eqaute the integrands, to obtain
$$\tag{3}\frac{\partial (\rho \mathbf{u})}{\partial t} +\nabla \cdot(\rho \mathbf{u} \mathbf{u}) =- \nabla p + \nabla \cdot \sigma $$
Now it is just a matter of simplifying both sides.
For the left-hand side we have
$$\tag{4}\frac{\partial (\rho \mathbf{u})}{\partial t} +\nabla \cdot(\rho \mathbf{u} \mathbf{u}) = \rho \left[\frac{\partial \mathbf{u}}{\partial t}+ \mathbf{u} \cdot \nabla \mathbf{u}\right] + \left[\frac{\partial \rho}{\partial t}+ \nabla \cdot (\rho\mathbf{u}) \right]\mathbf{u}$$
Conservation of mass is expressed in terms of the equation of continuity
$$\frac{\partial \rho}{\partial t}+ \nabla \cdot (\rho\mathbf{u}) =0,$$
which can be derived in a similar way. Substituting into (4) , we get
$$ \frac{\partial (\rho \mathbf{u})}{\partial t} +\nabla \cdot(\rho \mathbf{u} \mathbf{u}) = \rho \left[\frac{\partial \mathbf{u}}{\partial t}+ \mathbf{u} \cdot \nabla \mathbf{u}\right]= \rho \frac{D\mathbf{u}}{Dt}, $$
and substituting into (3) with this result,
$$ \rho \frac{D\mathbf{u}}{Dt} = -\nabla p + \nabla \cdot \sigma$$
It remains for you to expand $\nabla \cdot \sigma$ using the expression you gave for $\sigma$ in terms of the deformation tensor which should be $D = \nabla \mathbf{u} + [\nabla \mathbf{u}]^T.$