"but how can you have a "basis" that changes at every point?"
This is really the root of your problem. Mathematically (and physically) speaking, such a basis works fine. You just have a (hopefully temporary) conceptual problem.
Maybe try thinking of it this way: How does $\hat{x}$ know to point "to the right" in cartesian coordinates, at an arbitrary point in the coordinate system? You could define it as "parallel to the $x$-axis". Similarly for $\hat{y}$. Now how would you define $\hat{r}$? Parallel to a line joining the current point and the origin. $\hat{\theta}$? Perpendicular to $\hat{r}$. Not that different, right?
For the gradient, your mistake is that the components of the gradient vary contravariantly. On top of that, there is a issue with normalisation that I discuss below. I don't know if you are familiar with differential geometry and how it works, but basically, when we write a vector as $v^\mu$ we really are writing its components with respect to a basis.
In differential geometry, vectors are entities which act on functions $f : M \rightarrow \mathbb{R}$ defined on the manifold. Tell me if you want me to elaborate, but this implies that the basis vectors given by some set of coordinates are $\partial_\mu = \frac{\partial}{\partial x^\mu}$ and vary covariantly. Let's name those basis vectors $e_\mu$ to go back to the "familiar" linear algebra notation.
Knowing that, any vector is an invariant which can be written as $\vec{V} = V^\mu \partial_\mu$. The key here is that it is invariant, so it will be the same no matter which coordinate basis you choose.
Now, the gradient is defined in Euclidean space simply as the vector with coordinates $\partial_i f = \partial^if$ where $i = \{x,y,z\}$. Note that in cartesian coordinates covariant and contravariant components are the same. So, the invariant quantity is $\vec{\nabla}f = \partial^ife_i$. Note that, from what we did before, the components of a vector are to be treated as contravariant.
Now, since this expression is invariant we get, in general coordinates $\vec{\nabla}f = \partial^\mu fe_\mu$. So what you are looking for when computing the components is $\partial^\mu f = g^{\mu\nu}\frac{\partial f}{\partial x^\nu}$. This gives $\vec{\nabla}f = \frac{\partial f}{\partial r} e_r + \frac{1}{r^2}\frac{\partial f}{\partial \theta} e_\theta +\frac{1}{r^2 \sin^2\theta}\frac{\partial f}{\partial \phi} e_\phi$. This is still not what we're looking for.
This is due to the fact that the basis vectors are not normalised. Indeed, take a specific vector $e_I$. His components are $\delta^\mu_I$ by definition (it's a basis vector). Then, $|e_I|^2 = e_I^\mu e_I^\nu g_{\mu\nu} = g_{II}$. Great then ! Using $e_i' = e_i/\sqrt{g_{ii}}$ as normalised base vectors, we get the right answer : $\vec{\nabla}f = \frac{\partial f}{\partial r} e_r' + \frac{1}{r}\frac{\partial f}{\partial \theta} e_\theta' +\frac{1}{r \sin\theta}\frac{\partial f}{\partial \phi} e_\phi'$
Let's move onto the divergence. It seems easier since it's a scalar, there's no basis vector to mess around with. Well, actually there are still some problem attached to that. Your formula is right, again, except that when you write the invariant formula $\nabla_\mu V^{\mu}$ you implicitly use the basis that we defined earlier. This means that you are not working on a normalised basis. We know that the $\vec{V}$ you used is $\vec{V} = V^\mu e_\mu = V^\mu \sqrt{g_{\mu\mu}}e_\mu'$. So to compare formula you have to introduce the vector with respect to the normalised coordinates, $A^\mu= V^\mu\sqrt{g_{\mu\mu}}$. I'll let you plug it in your (correct) formula to see that it works.
To conclude, your formula for the curl should be right. Just be careful to use the right normalisations for the vectors and you should be fine (also be careful of the tensorial form of the levi-civita tensor, which involves the determinant of the metric). I have not the faith to do the calculations for you, but you definitely should try it to ensure yourself you understood it well.
P.S: Just for completeness, for the divergence there is a quite useful formula which is also used in Sean Caroll book : $\vec{\nabla}\cdot\vec{V} = \frac{1}{\sqrt{g}}\partial_\mu(\sqrt{g}V^\mu)$, useful when you're too lazy to compute Christoffels.
Best Answer
You were pretty close already. There is a handy table on Wikipedia for a variety of coordinate systems. But for the polar system:
$$ \vec{\nabla} \cdot \vec{U} = \frac{\partial U_r}{\partial r} + \frac{1}{r} \frac{\partial U_\theta}{\partial \theta} $$
and you can look up the curl in the same table.
These can be derived from the Cartesian definitions by considering the total differentials:
$$ dr = \frac{\partial r}{\partial x} dx + \frac{\partial r}{\partial y} dy$$
and
$$ d\theta = \frac{\partial \theta}{\partial x} dx + \frac{\partial \theta}{\partial y} dy$$
Or to put it another way, you end up with:
$$ \frac{\partial \phi}{\partial r} = \frac{\partial \phi}{\partial x}\frac{\partial x}{\partial r} + \frac{\partial \phi}{\partial y}\frac{\partial y}{\partial r}$$
which is just the chain-rule for derivatives. $\phi$ is any variable. A similar expression can be written for the $\theta$ coordinate. You can use this to go between any coordinate systems you would like to go between and relate it all back to the Cartesian system by successive chaining.
As for a physical insight, that's a little bit harder. But the way I like to look at it for the polar (aka cylindrical) system is that the unit length of $\theta$ changes as $r$ changes. If you draw a big circle vs a small circle, the arc length of the wedge defined by an angle $\theta$ is certainly different. And so it is natural that the derivative in the $\theta$ direction will have some radial dependence.