Your last suggestion is correct : the use which you are making of dimensional analysis is not justified.
Dimensional Analysis has two purposes : (1) to check that equations or terms in equations are commensurate; and (2) to find combinations of quantities with particular dimensions or no dimensions at all. It is not able to derive physically meaningful formulas, and certainly not a unique formula, which is what you seem to be expecting it to do.
The 2nd purpose can be achieved using Rayleigh's Method, or Buckingham's Pi Theorem which is a more formal version of it. For this purpose the method works, and it is easy to see why it works. That to me is sufficient justification.
In your 1st example, Rayleigh's Method tells you that the formula $\Delta P=k\rho gh$ (where $k$ is a dimensionless constant) is dimensionally consistent, but it cannot tell you that this formula is physically meaningful or correct for the application you have in mind. It cannot tell you what value of $k$ you should use. It can only tell you what dimensions $k$ must have for the formula $\Delta P=k\rho gh$ to be dimensionally consistent.
In the 2nd example, the method cannot tell you that the force of gravitational attraction between two masses should be $F=Gm_1m_2/r^2$. The method does not generate any dimensionless combinations of the quantities $F, m_1, m_2, r$. This is not a failure of Rayleigh's Method because no dimensionless combinations of these variables can be derived by any method. As with your 1st example, dimensional analysis can only tell you what dimensions the constant $G$ must have for the formula $F=Gm_1m_2/r^2$ to be dimensionally correct.
The strengths and limitations of the method are stated in the wikipedia article :
The Buckingham π theorem provides a method for computing sets of dimensionless parameters from given variables, even if the form of the equation remains unknown. However, the choice of dimensionless parameters is not unique; Buckingham's theorem only provides a way of generating sets of dimensionless parameters and does not indicate the most "physically meaningful".
Two systems for which the [dimensionless] parameter [are equal] are called similar. As with similar triangles, they differ only in scale. They are equivalent for the purposes of the [unknown] equation. The experimentalist who wants to determine the form of the equation can choose the most convenient [system to investigate].
Most importantly, Buckingham's theorem describes the relation between the number of variables and [the number of] fundamental dimensions.
Your intuition is wrong on this. Consider one-dimensional steady flow, say in the $x$-direction, with a velocity gradient in the $y$-direction. Thus the particles at a given level have average velocity $$\bar{\mathbf u}=(u(y), 0, 0)^T,$$ and fluctuating velocities $${\mathbf u}'=(u', v', w')^T.$$
Let's consider particles that at time $t_0$ are located at $(x,y_0)^T$, which have velocities ${\mathbf u}=(u_0+u', v', w')^T.$ These particles will, on average, travel a distance of the mean free path length $l$ at that velocity, before hitting other particles. The particles will thus have migrated to a different $y$ position, where the average particle velocity will be $$\bar{\mathbf u(y)}=(u(y_0)+(y-y_0)\frac{\partial u}{\partial y}, 0, 0)^T.$$ Notice that the average difference in the $x$-component of the velocity of such particles will therefore be proportional to the mean free path $l$ times an integral $I$ over the distribution of $v'$ and $w'$ velocities which does not matter here: We have $y-y_0=I\,l$. The mean velocity difference for such particles is therefore just $\bar{\Delta u}=I\,l\,(\partial u/\partial y)$.
Since the mean velocity is assumed to stay constant, such particles will have their velocity adjusted to the one at their new $y$-position. Viscous forces correspond to the work required to achieve this. These forces must therefore be proportional to the velocity gradient times the mean free path length.
P.S.: Also see the derivation in the Wikipedia article on viscosity.
Best Answer
A simple example: $x=at^2$ is dimensionally homogeneous, but the true equation is $x=\frac{1}{2}at^{2}\;\;$ (in the case $\;v_{0}=0,x_{0}=0$).