Let us clear up a few things first:
1. Bernoulli's equation is first and foremost a "first integral" of the equations of motion along a "streamline".
A streamline at a given instant is a curve in space whose tangent at a particular point equals the field velocity at that point (and instant).
2. In the "steady" flow, since the velocity field is locally independent of the time variable, the streamlines will coincide with particle trajectories, and since integrating the particle's equations of motion along its trajectory over a time interval gives us the energy balance, the standard Bernoulli equation does, indeed, COINCIDE with energy conservation.
3. Note, however, that the underlying DEFINITIONS are different:
For Bernoulli(-like) quantities, we integrate SPATIALLY, whereas for energy(-like) quantities, we integrate TEMPORALLY.
For unsteady flows, then, a "Bernoulli" procedure does NOT give an energy conservation picture, but rather
A picture of the INSTANTANTENOUS spatial distribution of dynamical quantities like pressure along curves that do NOT represent the trajectories of any particles whatsoever.
(Energy is not, in general, conserved along a strip of the streamline, the local time derivative term tells us how much is lost or gained on it)
To develop the maths sufficient to derive it, let \gamma be a variable along some streamline, and let the family of streamlines in the field at particular instant t be given as:
\vec{S}(\gamma;t,\vec{x}_{0}),\vec{S}(0;t,\vec{x}_{0})=\equiv\vec{x}_{0}
That the streamline tangent is everywhere parallell to the velocity field vector \vec{v}(\vec{x},t) is given by the defining equation:
\frac{\partial\vec{S}}{\partial\gamma}=\vec{v}{\vec{S},t)
Varying \gamma in \vec{S} moves us therefore along the same streamline, varying \vec{x}_{0} shifts us to a neighbouring streamline at the same instant t, whereas the quantity \frac{\partial\vec{S}}{\partial{t}} tells us how the streamlines changes over time.
Now, for our general derivation in the inviscid case, using constant density for convenience, we introduce index notation for clarity.
We call the streamline x_{i}(\gamma,t), i=1,2,3, and our equations of motion reads:
v_{i,t}+v_{j}v_{i,j}=-\frac{1}{\rho}p_{,i}-g_{i},
where indices after comma represent the variable we differentiate with respect to, and g_1=g_2=0.
Now, we take the dot product with \frac{\partial{x}_{i}}{\partial\gamma}=x_{i,\gamma}(=v_{i})
and gain (by making a legit switch of variable name):
v_{i,t}x_{i,\gamma}+(\frac{1}{2}v_{j}^{2})_{i}x_{i,\gamma}+\frac{1}{\rho}p_{,i}x_{i,\gamma}+gx_{3,\gamma}=0
Integrating the result between, say, \gamma_{0} to \gamma_{1} yields the already noted result in a previous post.
The remaining integral of the time-derivative of the velocity is usually difficult to estimate.
