I might have a proof, but I'm not sure if it's valid. I tried writing Phi as an unknown function or x, y and z, and u as a column vector of 3 more unknown 3D functions. Then I expanded Div[(Phi)(u)] using the product rule on each term - so I had a new column vector with two terms in each direction. Is that allowed, since each entry was a function, but not a vector?
Then I rearranged things to show that the result was equivalent to the proof desired, the result being applicable to n-dimensional space since since the operations involved were all independent of number of dimensions.
Does that work? It feels like I'm assuming things I shouldn't.
In any case, just to check, the later formula in the question, is correctly interpreted as stating that the rate of change of fluid density in a given volume is equal to the rate of flow of fluid in or out of the volume, no?