Let$a^2+b^2+c^2>R^2$, calculate $\int_{x^2+y^2+z^2\leq R^2}\frac{dxdydz}{\sqrt{(x-a)^2+(y-b)^2+(z-c)^2}}$.
Let $(x,y,z)=(a+r\sin\phi\cos \theta,b+r\sin\phi\sin \theta,z=c+r\cos\phi)$, but I don’t know how to transform the integration region. Is there any good idea?
Consider a transformation to a coordinate system centered at the point $(a,b,c)$ but rotated such that the origin lies on the $z$ axis relative to the point $(a,b,c)$. This gives us the integral
$$I = \iiint\limits_{x^2+y^2+\left(z-\sqrt{a^2+b^2+c^2}\right)^2\leq R^2}\frac{dV}{\sqrt{x^2+y^2+z^2}}$$
Then convert to spherical coordinates and do the angular integral first.
$$I = \int_0^{2\pi}\int_{\sqrt{a^2+b^2+c^2}-R}^{\sqrt{a^2+b^2+c^2}+R}\int_0^{\cos^{-1}\left(\frac{r^2+a^2+b^2+c^2-R^2}{2r\sqrt{a^2+b^2+c^2}}\right)}r\sin\theta\:d\theta\:dr\:d\theta$$
$$= \frac{\pi}{\sqrt{a^2+b^2+c^2}}\int_{\sqrt{a^2+b^2+c^2}-R}^{\sqrt{a^2+b^2+c^2}+R}2r\sqrt{a^2+b^2+c^2}-(r^2+a^2+b^2+c^2-R^2)\:dr$$
$$= \frac{\pi}{\sqrt{a^2+b^2+c^2}}\left(2R(a^2+b^2+c^2+R^2)-2R(a^2+b^2+c^2)-\frac{2}{3}R^3\right)$$
$$\implies \boxed{I = \frac{4\pi R^3}{3\sqrt{a^2+b^2+c^2}}}$$
which is exactly the volume of the sphere times the value of the original integrand at the origin.