In a curvilinear system with coordinates: $u_1, u_1, u_3$.
I want to understand why: $$\nabla\cdot \vec{V} \neq \frac{1}{h_1} \frac{\partial V_1}{\partial u_1} +\frac{1}{h_2} \frac{\partial V_2}{\partial u_2} + \frac{1}{h_3}\frac{\partial V_3}{\partial u_3}. \tag1$$
My reasoning to explain why eq.(1) is not true is the following:
Taking the first term of the dot product:
$$\frac{\hat{e_1} }{h_1}\frac{\partial}{\partial u_1} (V_1 \hat{e_1}) = \frac{\hat{e_1}}{h_1} \big( \hat{e_1} \frac{\partial V_1}{\partial u_1} + V_1 \frac{\partial \hat{e_1}}{\partial u_1} \big)$$
This seems rare to me since the definition of the dot product is:
$$A \cdot B = A_1B_1 + A_2B_2 + A_3B_3 $$
And using that definition the divergence should be eq.(1) (with an equality).
I would appreciate any references that explain this, since all I can find is the formula for the divergence in curvilinear systems but without a deduction.
Any good book on vector or tensor analysis should provide you with a first principles derivation for some specific curvilinear coordinate systems (spherical and cylindrical systems for example). For the more general case, you need to understand that your definition for the dot product of two vectors is valid only for Cartesian coordinates. Why is that? Because in general the dot product is defined as the multiplicative operation that converts two vector quantities into a scalar quantity. Likewise the divergence is defined as the first order differential operator that yields a scalar quantity when operating on a vector. In tensor analysis this operator is called the covariant derivative and for curvilinear coordinate systems the usual partial derivatives are augmented by terms called Christoffel symbols. These extra terms render the divergence a scalar quantity. This answer should provide you with keywords that will take you to a complete understanding.