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 (*).
Take the 2D, incompressible continuity equation
$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$$
and multiply it with $u$, leading to
$$u\frac{\partial u}{\partial x}+u\frac{\partial v}{\partial y}=0\,.$$
Now add the above equation to the reduced x-direction equation (the one you wrote correctly in your post); this will lead to
$$2u\frac{\partial u}{\partial x}+u\frac{\partial v}{\partial y}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho}\frac{dp}{dx}+\nu\frac{\partial^2 u}{\partial y^2}$$
$$\Rightarrow \frac{\partial u^2}{\partial x}+\frac{\partial(uv)}{\partial y}=-\frac{1}{\rho}\frac{dp}{dx}+\nu\frac{\partial^2 u}{\partial y^2} \,.$$
Notice that the pressure derivative is a total one, since there is no dependence on other variables.
I think this is the equation you're looking for, check if you can that no typo was made while copying the equation from Lienhard's textbook. Let me know if I missed anything!
Best Answer
You are conflating incompressible with inviscid. Here we have flow of an incompressible viscous fluid and the term $\nu \nabla^2 \mathbf{v}$ may not be neglected.
You are asked to find "a solution" under the assumptions that the flow is steady and uniform in the x-direction. This implies fully-developed unidirectional flow where the only non-vanishing component of velocity is $u$ (i.e., the component in the x-direction).
For an incompressible fluid, the equation of continuity gives
$$\nabla \cdot \mathbf{v} = \frac{\partial u}{\partial x} = 0$$
Thus, $u$ is a function only of the y-coordinate and the Navier-Stokes equations reduce to
$$0 = -\frac{1}{\rho}\frac{\partial p}{\partial x} +g \sin \alpha + \nu \frac{d^2u}{dy^2}, \\0 = -\frac{1}{\rho}\frac{\partial p}{\partial y} -g \cos \alpha $$
Since the flow is gravity-driven, we can neglect the x-component of the pressure pressure gradient $\frac{\partial p}{\partial x}$ in the first equation and find $u$ by applying the two boundary conditions to the solution of the second-order differential equation
$$\nu \frac{d^2 u}{dy^2} = -g \sin \alpha$$
The second equation allows for solution of the pressure as a function of the y-coordinate.
See if you can finish now.