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 Navier-Stokes equations apply as the fluid is Newtonian and, hence, viscous.
Neglecting body force terms other than gravity which can be absorbed into the pressure field, the $x$-component of velocity satisfies
$$\tag{*}\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z} = - \frac{1}{\rho} \frac{\partial p}{\partial x} + \frac{\mu}{\rho}\nabla^2 u,$$
where $\mu$ is the viscosity.
The flow is steady, fully developed ($u$ depends only on $y$), and unidirectional. Consequently the partial derivatives of $u$ with respect to $t$, $x$, and $z$ vanish and $v = w = 0$. Thus, (*) reduces to
$$\tag{**} \frac{d ^2u}{dy^2} = \frac{1}{\mu}\frac{\partial p}{\partial x} = \frac{G}{\mu}$$
The pressure gradient $G$ is independent of $y$. This follows from the $y$-component of the Navier-Stokes equations which reduces to
$$\frac{\partial p}{\partial y} = 0$$
Integrating both sides of (**) twice with respect to $y$ we get the quadratic function
$$u(y) = \frac{G}{2\mu}y^2 + C_1y + C_2,$$
wher $C_1$ and $C_2$ are integration constants.
Best Answer
Aside from an algebraic error, you have solved for the critical value of sound speed where the fluid attains Mach 1 -- hence, you have done nothing wrong. The sound speed will assume different values throughout the flow as a function of local temperature. Nevertheless the fluid will attain a condition where the velocity equals a specific sound speed (Mach 1) that depends on the density where the fluid is at rest. There are other constraints to isentropic flow, for example, a maximum attainable velocity which corresponds to the pressure falling to $0$.
At Mach 1 the fluid velocity equals the critical speed of sound $c_*$at local conditions, so $|\nabla\phi|^2 = c_*^2$ and substituting in the Bernoulli equation we have
$$\frac{c^2}{2} + \frac{k\gamma \rho^{\gamma-1}}{\gamma-1} = \frac{c_*^2}{2} + \frac{c_*^2}{\gamma-1} = C \\ \implies c_*^2 \frac{\gamma +1 }{2(\gamma -1)}= C \\ \implies c_*^2 = 2C \frac{\gamma-1}{\gamma+1} \quad (*)$$.
The constant C can be related to the speed of sound $c_0$ when the fluid is at rest. Under such stagnation conditions the fluid velocity is $0$ and, again, using the Bernoulli equation we have
$$ \frac{k\gamma \rho_0^{\gamma-1}}{\gamma-1} = \frac{c_0^2}{\gamma-1} = C $$
Substituting for C in (*) we obtain
$$c_*^2 = 2\frac{c_0^2}{\gamma-1} \frac{\gamma-1}{\gamma+1} = \frac{2c_0^2}{\gamma+1}, $$
which yields
$$\frac{c_*}{c_0} = \sqrt{\frac{2}{\gamma+1}}$$
Typical values for air are $\gamma = 1.4$ and $c_* \approx 0.913 c_0$.