Poiseuille flow: can I avoid the simplification that $u=u(y)$? Let's consider a viscous, incompressible and irrotational fluid that flows in a canal in two dimensions. It is reasonable to assume that the velocity of the fluid has the form $\vec v = u(x,y,t) \hat x $, or in other words, we are assuming that the velocity has just the horizontal component $u(x,y,t)$.
Under these hypotheses, the incompressibility condition $\vec \nabla \cdot \vec v = 0$ gives us that $\partial_x u = 0$, i.e. $u$ depends only on $t$ and $y$.
The Navier-Stokes equations then become
$$ (1) \qquad \begin{cases}
\partial_tu = -\partial_xp + \nu \partial_{yy}u \\
-\partial_y p = 0
\end{cases}$$
At this point one usually says "suppose that the velocity does not depend on time": the first equation of the system $(1)$ simplifies to $\nu \partial_{yy}u = \partial_x p $ and one obtain a closed form of the pressure $p$ and the velocity $u$.
My question is: while the assumptions made at the beginning seem resonable with the usual context of fluid dynamics, the last hypothesis seems to be made just to derive in few steps a solution of the problem. Although it is said "a physical theory will always be wrong" because of the assumptions made to translate the real world into some mathematical equations, I ask:

Can the system $(1)$ have solutions without the simplification that $u = u(y)$?

 A: You have reduced the Navier-Stokes equations to a simple form by making a number of assumptions that may or may not be realistic or commensurate with a physical configuration.  To begin with, restricting to two-dimensional flow where there is no dependence on the spanwise coordinate $z$ nor a non-zero velocity component $w$ in that direction can only be true assuming infinite extent of the domain in the $z-$direction.  Furthermore, by assuming there is no dependence on the streamwise coordinate $x$, you are making the assumption of fully developed flow far from an entrance to the channel where the velocity field must assume a more complicated form. Nevertheless, if the physical dimensions are such that the the spanwise extent is very large and we are examining the flow far from an entrance, this can be a good approximation.
Of more immediate concern is that you have not mentioned boundary conditions.  If the channel is confined between solid stationary boundaries at $y = 0,L$, then we must impose the boundary conditions $u(0) = u(L) =0$.  A further assumption is needed about the pressure gradient.  The simplest would be that it is constant, $\frac{\partial p}{\partial x} = -G < 0$. In this case, (1) admits a steady solution obtained by first integrating twice with respect to $y$ to get
$$u(y) = C_1 +C_2y - \frac{Gy^2}{2\nu},$$
and, after applying boundary conditions, we obtain the velocity profile for steady plane Poiseuille flow:
$$u(y) = \frac{GL}{2\nu}\left( y - \frac{y^2}{L}\right)$$
Again, this is an idealization for the flow between infinite parallel plane boundaries with a constant pressure gradient.
Changing the configuration or boundary conditions by introducing time-dependency will lead to solutions of (1) with the full functional form $u(y,t)$.  An example would be when the pressure gradient is itself time-dependent, for example oscillating sinusoidallly in time.
A classic variation with a time-dependent solution is the Stokes problem. In this case the domain is of semi-infinite extent $0 \leqslant y < \infty$ and the flow is driven not by pressure but, instead, by the sinusoidal oscillation of a plate at $y=0$ requiring the boundary condition $u(0,t) = U\cos \omega t$.  In this case, the pressure gradient is taken to be zero and the closed form solution to (1) is
$$u(y,t) = Ue^{-\sqrt{\frac{\omega}{2\nu}}y}\cos \left(\omega t - \sqrt{\frac{\omega}{2\nu}}y \right)$$
A closed-form time dependent solution can even be obtained when $0 \leqslant y \leqslant L$ and there is a stationary wall at $L$ where $u(L,t) =0$ for all $t$.
Even in simple configurations like flow between infinite parallel plates,  the full Navier-Stokes equations admit extremely complex  time-dependent solutions that reflect the instability of unidirectional, laminar flow and the transition to fully developed turbulence.
