I'm assuming that you already know how to get the curl for a vector field in Cartesian coordinate system. When you try to derive the same for a curvilinear coordinate system (cylindrical, in your case), you encounter problems. Cartesian coordinate system is "global" in a sense i.e the unit vectors $\mathbb {e_x}, \mathbb {e_y}, \mathbb {e_z}$ point in the same direction irrepective of the coordinates $(x,y,z)$. On the other hand, the curvilinear coordinate systems are in a sense "local" i.e the direction of the unit vectors change with the location of the coordinates.
For example, in a cylindrical coordinate system, you know that one of the unit vectors is along the direction of the radius vector. The radius vector can have different orientation depending on where you are located in space. Hence the unit vector for point A differs from those of point B, in general.
I'll first try to explain how to go from a cartesian system to a curvilinear system and then just apply the relevant results for the cylindrical system. Let us take the coordinates in the new system as a function of the original coordinates. $$q_1 = q_1(x,y,z) \qquad q_2 = q_2(x,y,z) \qquad q_3 = q_3(x,y,z) $$
Let us consider the length of a small element
$$ ds^2 = d\mathbf{r}.d\mathbf{r} = dx^2 + dy^2 + dz^2 $$
The small element $dx$ can be written as
$$ dx = \frac{\partial x}{\partial q_1}dq_1 + \frac{\partial x}{\partial q_2}dq_2+\frac{\partial x}{\partial q_3}dq_3 $$
Doing the same for $dy$ and $dz$, we can get the distance element in terms of partial derivatives of $x,y,z$ in the new coordinate system. This will be of the form
$$ ds^2 = \sum_{i,j} \frac{\partial \mathbf{r}}{\partial q_i} .
\frac{\partial \mathbf{r}}{\partial q_j} dq_i dq_j = \sum_{i,j} g_{ij} dq_i dq_j $$
Here $\frac{\partial \mathbf{r}}{\partial q_j}$ represents the tangent vectors for $q_i = $constant , $i\neq j$. For an orthoganal coordinate system, where the surfaces are mutually perpendicular, the dot product becomes
$$ \frac{\partial \mathbf{r}}{\partial q_i} .
\frac{\partial \mathbf{r}}{\partial q_j} = c \delta_{ij}$$
where the scaling factor $c$ arises as we haven't considered unit vectors. These factors are taken as $$c = \frac{\partial \mathbf{r}}{\partial q_i} .
\frac{\partial \mathbf{r}}{\partial q_i} = h_i^2$$
$$ ds^2 = \sum_{i=1}^3 (h_i dq_i)^2 = \sum_i ds_i^2$$
Hence the length element along direction $q_i$ is given by $ds_i = h_i dq_i $.
Now, we are equipped to get the curl for the curvilinear system. Consider an infinitesimal enclosed path in the $q_1 q_2$ plane. And evaluate the path integral of the vector field $\mathbf{V}$ along this path.
$$\oint \mathbf{V}(q_1,q_2,q_3).d\mathbf{r} = \oint \mathbf{V}.\left( \sum_{i=2}^3 \frac{\partial \mathbf{r}}{\partial q_i} dq_i\right)$$
(q1, q2+ds_2) (q1+ds_1, q2+ds_2)
-----------<------------
| |
| |
V ^
| |
|---------->-----------|
(q1, q2) (q1+ds_1, q2)
$$ \oint \mathbf{V}.d\mathbf{r} = V_1 h_1 dq_1 -
\left( V_1 h_1 + \frac{\partial V_1 h_1}{\partial q_2} dq_2\right)dq_1
- V_2 h_2 dq_2 +
\left( V_2 h_2 + \frac{\partial V_2 h_2}{\partial q_1} dq_1\right)dq_2$$
$$ = \left( \frac{\partial V_2 h_2}{\partial q_1} -
\frac{\partial V_1 h_1}{\partial q_2} \right)dq_1 dq_2$$
From Stokes theorem,
$$ \oint \mathbf{V}.d\mathbf{r} =
\int_S \nabla \times \mathbf{V} . d\mathbf{\sigma} =
\nabla \times \mathbf{V} . \mathbf{\hat{q_3}} (h_1 dq_1) (h_2 dq_2)
= \left( \frac{\partial V_2 h_2}{\partial q_1} -
\frac{\partial V_1 h_1}{\partial q_2} \right)dq_1 dq_2 $$
Hence the $3$ component of the curl can be written as
$$(\nabla \times \mathbf{V})_3 = \frac{1}{h_1 h_2} \left( \frac{\partial V_2 h_2}{\partial q_1} -\frac{\partial V_1 h_1}{\partial q_2} \right) $$
Similarly, other components can be evaluated and all the components can be assembled in the familiar determinant format.
$$\nabla \times \mathbf{V} = \frac{1}{h_1 h_2 h_3} \begin{vmatrix}
\mathbf{\hat{q_1}}h_1 & \mathbf{\hat{q_2}}h_2 & \mathbf{\hat{q_3}}h_3\\
\frac{\partial}{\partial q_1} & \frac{\partial}{\partial q_2} & \frac{\partial}{\partial q_3} \\
V_1 h_1 & V_2 h_2 & V_3 h_3 \\
\end{vmatrix}$$
Now the expression for the curl is ready. All we need to do is find the values of $h$ for the cylindrical coordinate system. This can be obtained, if we know the transformation between cartesian and cylindrical polar coordinates.
$$ (x,y,z) = (r\cos\phi, r\sin\phi, z)$$
Now the length element
$$ ds^2 = dx^2 + dy^2 + dz^2 = (d(r\cos\phi))^2 + (d(r\sin\phi))^2 + dz^2 $$
Simplifying the above expression, we get
$$ ds^2 = (dr)^2 + r^2(d\phi)^2 + (dz)^2 $$
From the above equation, we can obtain the scaling factors, $h_1 = 1$ , $h_2 = r$, $h_3 = 1$. Hence the curl of a vector field can be written as,
$$ \nabla \times \mathbf{V} = \frac{1}{r}
\begin{vmatrix}
\mathbf{\hat{r}} & r\mathbf{\hat{\phi}}& \mathbf{\hat{z}}\\
\frac{\partial}{\partial r} & \frac{\partial}{\partial \phi} & \frac{\partial}{\partial z} \\
V_r & r V_\phi & V_z \\
\end{vmatrix}
$$
"...I feel like this is too good to be true". Your question is perfectly legit since changes of coordinates are always tricky for differential operators. But luckily there is a general and well known theory, dubbed "curvilinear coordinates".
The answer is: YES, it is correct, but for a scalar field only. Do not try to extend naively this formula to higher-rank fields. Example: the directional derivative of a vector field in cylindrical coordinates is much more complex than this. Take a look at:
https://link.springer.com/content/pdf/bbm%3A978-3-0348-8579-9%2F1.pdf
The general topic of your question is how "orthogonal curvilinear coordinates" work: this includes the standard cartesian, cylindrical and spherical coordinates.
Best Answer
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). $$