Consider a small volume element
$$\{(s_1,s_2,s_3)\in \mathbb{R}^3: x_1 \leq s_1 \leq x_1 + \Delta x_1,x_2 \leq s_2 \leq x_2 + \Delta x_2,x_3 \leq s_3 \leq x_3 + \Delta x_3\}.$$
The surface (stress) force acting on that element in the $i$ direction consists of contributions from each of three pairs of faces where coordinates $s_1$, $s_2$ and $s_3$ are constant, respectively.
Consider one pair of where $s_1$ is fixed: $s_1 = x_1$ and $s_1 = x_1 + \Delta x_1$. The net force in the $i$ direction is
$$F_{i,1}\approx \sum_{j=1}^{3}\sigma_{ij}n_j|_{s_1=x_1 + \Delta x_1}\Delta x_2\Delta x_3+\sum_{j=1}^{3}\sigma_{ij}n_j|_{s_1=x_1}\Delta x_2\Delta x_3$$
This is approximate because we evaluate $\sigma_{ij}(s_1,s_2,s_3)$ at representative coordinates $s_2= x_2$ and $s_3=x_3$ . In general, we assume that the stress tensor is a continuously differentiable function of spatial position.
Here we use the fact that the surface force (in terms of the stress tensor and normal vector) is the product $\mathbb{\sigma}\cdot\mathbb{n}.$ On the face $s_1 = x_1 + \Delta x_1$ we have $\mathbb{n} = (1,0,0)$, and on the face $s_1 = x_1$ we have $\mathbb{n} = (-1,0,0).$
Hence,
$$F_{i,1}\approx \sigma_{i1}|_{s_1=x_1 + \Delta x_1}\Delta x_2\Delta x_3-\sigma_{i1}|_{s_1=x_1}\Delta x_2\Delta x_3$$
The contribution to the force per unit volume is
$$\frac{F_{i,1}}{\Delta x_1 \Delta x_2 \Delta x_3}\approx \frac{\sigma_{i1}|_{s_1=x_1 + \Delta x_1}-\sigma_{i1}|_{s_1=x_1}}{\Delta x_1}.$$
Now consider the limit as the element shrinks to the point $(x_1,x_2,x_3)$ and $\Delta x_k \rightarrow 0.$
Then the contribution to the net surface force per unit volume at that point is
$$\hat{F_{i,1}}= \lim_{\Delta x_1 \rightarrow 0}\frac{\sigma_{i1}|_{s_1=x_1 + \Delta x_1}-\sigma_{i1}|_{s_1=x_1}}{\Delta x_1}=\frac{\partial \sigma_{i1}}{\partial x_1}.$$
Applying the same derivation to the other pairs of faces we obtain the total contribution to the net surface force per unit volume (in the $i$ direction) as
$$\hat{F_{i}}=\sum_{j=1}^{3}\frac{\partial \sigma_{ij}}{\partial x_j}$$
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 (*).
Best Answer
You are on the right track, but it would be easier if you combined the second term into a single tensor (the rate-of-strain tensor) rather than treating it as multiple terms.
A good reference on this subject is Rutherford Aris's book titled "Vectors, Tensors, and the Basic Equations of Fluid Mechanics". It is available in paperback and is my go-to reference for things like this. Page 180 gives the general equation for the contravariant components of the stress tensor of a Newtonian fluid as
$$ \sigma^{ij} = ( - p + \lambda e^m{}_m ) g^{ij} + 2 \mu e^{ij} \,, $$
where the covariant components of the rate-of-strain tensor are defined as
$$ e_{ij} \equiv \tfrac{1}{2} ( u_{j,i} + u_{i,j} ) \,. $$
Here the commas just refer to the covariant derivative. It seems that in your case, the flow is incompressible, so $e^m{}_m = 0$, and the 2nd term drops out leaving
$$ \sigma^{ij} = - p g^{ij} + 2 \mu e^{ij} \,, $$
This basically is your answer but written in terms of the contravariant components and using the covariant derivative.
Unfortunately most books either use vector notation or use Cartesian index notation (like G. K. Batchelor's "Introduction to Fluid Mechanics" or Landau/Lifshitz's book). Aris's book is one of the few that I know that uses general tensor notation, so it is a good reference for this topic.