It's not hard to derive the equation for the Del operator in cylindrical coordinates from the Del operator in cartesian coordinates. From $$\nabla = \hat{\bf{x}} \frac{\partial}{\partial x} + \hat{\bf{y}}\frac{\partial}{\partial y} + \hat{\bf{z}}\frac{\partial}{\partial z}$$ we arrive at $$\nabla^* = \hat{\bf{r}} \frac{\partial}{\partial r} + \mathbf{\hat{\theta}} \frac{1}{r} \frac{\partial}{\partial \theta} + \hat{\bf{z}} \frac{\partial}{\partial z}$$ using standard relationships between $(x,y,z)$ and $(r,\theta,z)$.
What I'm wrestling with is how the dot product operation is defined with the $\nabla^*$ operator. That is for $\textbf{v} = <v_r, v_\theta, v_z>$ we define
\begin{equation}\nabla^* \cdot \textbf{v} = \frac{1}{r}\frac{\partial (r v_r)}{\partial r} + \frac{1}{r} \frac{\partial v_{\theta}}{\partial \theta}+ \frac{\partial v_z}{\partial z} \end{equation}
The latter terms are consistent with my understanding of the dot product in cartesian coordinates but the first term has factors of $\frac{1}{r}$ and $r$ that I don't immediately intuit.
Deriving the continuity equation in cylindrical coordinates for fluids of constant properties using a control volume argument also brings in these factors of $r$, arising from the fact that fluid leaving in the radial direction moves through a cross sectional area of $(r+dr)d\theta dz$ which is larger than what it enters through of area $rd\theta dz$.
I skimmed through Wolfram's definition at equations (112)-(115) and I suspect that the discrepancy arises from the Christoffel symbols (42)-(47) and more basically from the gradient components given at (34),(37).
Can anyone inspire some intuition for me here? Or is this a hopelessly "can't see the forest through the trees" situation?
I'm going to give a rough outline that you can use to remember this, rather than a derivation. In an orthogonal coordinate system, the scale factors determine all the formulae for differential operators: these come from the definition $$ \frac{\partial \mathbf{r}}{\partial q_i} = h_i \mathbf{e}_i, $$ where $\mathbf{e}_i$ is a unit vector, and $h_i$ is a function taken as positive. We set it up so that the $\mathbf{e}_i$ at a point are orthogonal (this is actually what it means to say the coordinate system is orthogonal), $$\mathbf{e}_i \cdot \mathbf{e}_i = \delta_{ij}.$$ Then the line element looks like $$ d\mathbf{r} = \sum_i \frac{\partial \mathbf{r}}{\partial q_i} \, dq_i = \sum_i \mathbf{e}_i h_i \, dq_i. $$ Then we can define the gradient by using the definition $df = \operatorname{grad}{f} \cdot d\mathbf{r}$ (think in Cartesians: this makes sense by the chain rule formula): $$ df = \sum_i \frac{\partial f}{\partial q_i} \, dq_i = \sum_i \frac{1}{h_i}\frac{\partial f}{\partial q_i} h_i \, dq_i = \sum_{i,j} \frac{1}{h_i}\frac{\partial f}{\partial q_i} \mathbf{e}_i \cdot \mathbf{e}_j h_j \, dq_j = \left( \sum_i \frac{\mathbf{e}_i}{h_i} \frac{\partial f}{\partial q_i} \right) \cdot d\mathbf{r}, $$ using the orthogonality of the $\mathbf{e}_i$, so $$ \operatorname{grad}f = \sum_i \frac{\mathbf{e}_i}{h_i} \frac{\partial f}{\partial q_i}. $$
Right, now the divergence. Let's think about what we want. My way of remembering what the divergence does is to think of it as a sort of formal adjoint to the gradient operator, because $$ \operatorname{div}(f\mathbf{u}) = \operatorname{grad}{f} \cdot \mathbf{u} + f \operatorname{div}{\mathbf{u}}, $$ so the divergence theorem gives $$ \int_V \operatorname{grad}{f} \cdot \mathbf{u} \, dV = - \int_V f \operatorname{div}{\mathbf{u}} \, dV + \text{(boundary terms)}. $$ The volume element is a cuboid spanned by three line elements, so $dV = h_1 h_2 h_3 \, dq_1 \, dq_2 \, dq_3 $.
Let's just look at one term: the others will work in the same way by the symmetry in the operators and linearity: suppose $u = u_1 \mathbf{e}_i$, and then $$ \int_V \operatorname{grad}{f} \cdot \mathbf{u} \, dV = \iiint \frac{1}{h_1}\frac{\partial f}{\partial q_i} \mathbf{e}_i \cdot \mathbf{e}_1 u_1 \, h_1 h_2 h_3 \, dq_1 \, dq_2 \, dq_3 \\ = \iiint \frac{\partial f}{\partial q_1} u_1 \, h_2 h_3 \, dq_1 \, dq_2 \, dq_3 $$ Now, apply integration by parts to the $q_1$ integral: $$ \iiint \frac{\partial f}{\partial q_1} u_1 \, h_2 h_3 \, dq_1 \, dq_2 \, dq_3 = \text{(boundary terms)} - \iiint f \frac{\partial}{\partial q_1} (h_2 h_3 u_1) \, dq_1 \, dq_2 \, dq_3 $$ Okay, but now I want to get back to a $dV$, so multiply and divide by $h_1 h_2 h_3$ to obtain $$ \iiint \frac{\partial f}{\partial q_1} u_1 \, h_2 h_3 \, dq_1 \, dq_2 \, dq_3 = \text{(boundary terms)} - \iiint f \left( \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial q_1} (h_2 h_3 u_1) \right) \, dV \\ = \text{(boundary terms)} - \iiint f \operatorname{div} \mathbf{u} \, dV, $$ so we conclude that the divergence is given by the expression $$ \operatorname{div} \mathbf{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) $$
Right, fine. So what is the divergence in cylindrical coordinates, like you actually asked? Well, we just have to find the scale factors. $\mathbf{r} = (r\cos{\theta},r\sin{\theta},z)$, so $$ h_r \mathbf{e}_r = 1(\cos{\theta},\sin{\theta},0)\\ h_{\theta} \mathbf{e}_{\theta} = r(-\sin{\theta},\cos{\theta},0)\\ h_z \mathbf{e}_z = 1(0,0,1), $$ so the only nontrivial $h$ is $h_{\theta}=r$. Thus $$ \operatorname{div} \mathbf{u} = \frac{1}{r} \left( \frac{\partial}{\partial r} (r u_r) + \frac{\partial}{\partial \theta} (u_\theta) + \frac{\partial}{\partial z} (r u_z) \right) = \frac{1}{r} \frac{\partial}{\partial r} (r u_r) + \frac{1}{r}\frac{\partial}{\partial \theta} (u_\theta) + \frac{\partial}{\partial z} (u_z). $$