If you really want to use the divergence theorem, the piece of the ball that you have to remove (and calculate the flux through) is almost like a cylinder of radius along the z-axis, with curved top and bottom. Not a sphere with radius .

Through the side of a cylinder of radius , the flux is easy to compute since the normal vector is in cylinder coordinates, so that you just need to integrate over the side: and, as tends to 0, and .

The problem is that the top and bottom sides remain, but you can prove that the flux through them tends to 0 as tends to 0. For that, it suffices to prove that can be continuously extended at the poles of the sphere (where ), because then the flux is asymptotically equivalent, as , to the value of this function at the pole times the area of the piece of sphere we consider, which converges to 0. To prove the continuous extension, you can show that where is the latitude (defined in the sphere); this is simple trigonometry if you draw an appropriate sketch.

...In fact, in such a symmetric situation, with such a nice vector field, and a singularity that prevents you to nicely apply the divergence theorem, it seems way simpler to me to directly use the definition of the flux as an integral over the surface.

By the above computation, the scalar product equals . And you need to integrate this on the sphere. The easiest way is to use spherical coordinates in the sphere: this gives where and as expressed above (depending on only). This is because is the usual surface element in spherical coordinates.