Expression for the Volume when Deriving Curvilinear Divergence

Question

I am trying to derive the expression for divergence of a vector field, $\vec{u}(\vec{r})$, in orthogonal curvilinear coordiantes $q_{1},q_{2},q_{3}$ with scale factors $h_{1},h_{2},h_{3}$ \begin{equation} \vec{\nabla}\cdot\vec{u}=\frac{1}{h_{1}h_{2}h_{3}}\left[\frac{\partial}{\partial q_{1}} (h_{2}h_{3}u_{1})+\frac{\partial}{\partial q_{2}}(h_{3}h_{1}u_{2})+\frac{\partial}{\partial q_{3}}(h_{1}h_{2}u_{3})\right] \end{equation} from the integral expression for divergence \begin{equation} \vec{\nabla}\cdot\vec{u}=\lim_{V\rightarrow0}\frac{1}{V}\oint_{S}\vec{u}\cdot d\vec{S} \end{equation} for a volume $V$ bounded by a surface $S$. I am broadly following this (from Equation 11). Doing algebra directly with objects such as $dq_{1}$ makes me uncomfortable so I am trying to follow the same outline but keeping everything expressely finite as Taylor expansions then at the end when the limit is taken higher order terms will give zero. As in the article I consider the flux throught the faces of a volume bounded by faces with one coordiante held constant. Then for each coordiate $q_{i}$ there are two faces each with that coordinate constant, one face with $q_{i}=a_{i}$ and another with $q_{i}=a_{i}+\delta q_{i}$. For the faces with the first coordinate, $q_{1}$, held constant I get a Taylor series whose first term is the expression the article has (but with $\delta$'s instead of $d$'s) \begin{equation} \int_{\text{surfaces}\perp q_{1}}\vec{u}\cdot d\vec{S}=\delta q_{1}\delta q_{2}\delta q_{3}\frac{\partial(h_{2}h_{3}u_{1})}{\partial q_{1}}+\text{higher order} \end{equation} where everything is evaluated at the point $(a_{1},a_{2},a_{3})$, and similarly for the other pairs of faces. I am trying to get an expression for the volume. The article simply divides by $dV=h_{1}h_{2}h_{3}dq_{1}dq_{2}dq_{3}$ but I am not convinced that is the right thing to divide by as that is the volume of a differential volume element at $q_{1},q_{2},q_{3}$ not some arbitrary volume (whose sides are surfaces of constant coordinates) and at this point we have not even taken the limit $V\rightarrow0$. I would like some sort of Taylor expansion for the volume such that the dominant term in the denominator when the limit of the ratio of flux and volume is taken is $V=h_{1}h_{2}h_{3}\delta q_{1}\delta q_{2}\delta q_{3}$ but I am not entirely sure how to go about this. I could write \begin{gather} V=h_{1}h_{2}h_{3}\delta q_{1}\delta q_{2}\delta q_{3}+\delta q_{1}\frac{\partial}{\partial q_{1}}h_{1}h_{2}h_{3}\delta q_{1}\delta q_{2}\delta q_{3}+\delta q_{2}\frac{\partial}{\partial q_{2}}h_{1}h_{2}h_{3}\delta q_{1}\delta q_{2}\delta q_{3}+\delta q_{3}\frac{\partial}{\partial q_{3}}h_{1}h_{2}h_{3}\delta q_{1}\delta q_{2}\delta q_{3}+\text{higher order} \end{gather} where everything is evaluated at the point $(a_{1},a_{2},a_{3})$. This 'feels' right (and ultimately leads to the right answer) but I can't explain in a clear way how one would arrive at that expression, or even if that is a correct expression.

Answer 1

The net flux across the two surfaces defined by $q_1=a, q_1=b$ whose faces are defined by a rectangle $R$ in the coordinates $q_2, q_3$ is exactly

Flux=$\int_{(q_2, q_3)\in R} \ [(u_1 h_2 h_3)|_{q_1=b} - (u_1 h_2 h_3)|_{q_1=a} ] dq_2 dq_3$

which by the one-variable FTC can be written as

Flux = $\int\int \int \ \frac{\partial}{\partial q_1} [u_1 h_2 h_3] dq_1 dq_2 dq_3$

Now make the substitution $dq_1 \ dq_2 \ dq_3= \frac{1}{h_1 h_2 h_3} dV$ and you will have an exact formula for the flux, written asa triple integral.

Do the same for the other pairs of faces and add the results together to get the total net flux expressed as a single triple integral of a (continuous) integrand. Then use the Mean Value Theorem for the triple integral to conclude that it equals the value of the integrand at some point within the small test region. As the test region tends to a point, this identity converges to the desired result.