Triple Integral $ \iiint\limits_S\frac1{\sqrt{\left(x-a\right)^{2}+\left(y-b\right)^{2}+\left(z-c\right)^{2}}}\;dx\;dy\;dz $ I used spherical coordinates and couldn't calculate this triple integral.
$$
\iiint_S\frac1{\sqrt{\left(x-a\right)^{2}+\left(y-b\right)^{2}+\left(z-c\right)^{2}}}\;dx\;dy\;dz
$$
where $S$ is a solid sphere of radius $R$
and center at the origin, and $\left ( a,b,c \right )$ is a fixed point outside this sphere.
 A: Method 1
This integral is actually very easy to evaluate because for any $\vec{p} \in \mathbb{R}^3$,
the function $\phi(\vec{x}) = \frac{1}{|\vec{x}-\vec{p}|}$ satisfies the 
Laplace's equation away from $\vec{p}$.
$$\nabla^2 \phi(\vec{x}) = 0 \quad\text{ for } x \in \mathbb{R}^3\setminus \{\vec{p}\}$$
This sort of function is called Harmonic function and one of its most useful property is the so called mean value property:

If $\phi$ is a function harmonic over a ball $B(\vec{q},R)$ centered
  at $\vec{q}$ with radius $R$, then
$$\phi(\vec{q}) = \frac{1}{4\pi R^2} \int_{\partial B(\vec{q},R)}
\phi(\vec{x}) dS(\vec{x}) = \frac{3}{4\pi R^3} \int_{B(\vec{q},R)}
\phi(\vec{x}) dV(\vec{x})$$

i.e. the average of $\phi(\vec{x})$ over the surface and the body of the ball $B(\vec{q},R)$
are both equal to the value of $\phi(\vec{x})$ at the center. Apply this to the integral
at hand, one can read off the value of the integral as:
$$\frac{4\pi R^3}{3\sqrt{a^2+b^2+c^2}}$$
Method 2
If one insists to evaluate the integral using polar-coordinates, one way to do it is to use
the so called multipole expansion.

Let  $\vec{p} = (a,b,c)$ and $\theta$ be the angle between the two
  vectors $\vec{x}$ and $\vec{p}$.
  Let $r_{<} = \min( |\vec{x}|, |\vec{p}| )$ and $r_{>} = \max(|\vec{x}|,|\vec{p}|)$, then
  $$\frac{1}{|\vec{x} - \vec{p}|} = \sum_{n=0}^{\infty}\frac{r_{<}^n}{r_{>}^{n+1}} P_n(\cos\theta)$$ 
  where $P_n(\cos\theta)$ are the Legendre's
  polynomials.

In our cases, $\vec{p}$ is outside the sphere, the triple integral can be rewritten as
$$\sum_{n=0}^{\infty} \frac{2\pi}{|\vec{p}|^{n+1}} \left( \int_0^R r^{n+2} dr \right) \left( \int_{0}^{\pi} P_n(\cos\theta) \sin\theta d\theta \right)$$
Using 
$$\int_{0}^{\pi} P_n(\cos\theta)\sin\theta d\theta = \begin{cases}2,&n = 0\\0,&n > 0\end{cases}$$
one can evaluate the triple integral as:
$$\frac{4\pi}{|\vec{p}|} \left( \int_0^R r^2 dr \right)
= \frac{4\pi R^3}{3\sqrt{a^2+b^2+c^2}}$$ 
Method 3
If one really want to evaluate the integral in the most elementary manner, then one
should use rotational symmetry to simplify the problem. We will assume $(a, b, c) = (0,0,-p)$ and in terms of cylindrical polar coordinate, the triple integral becomes:
$$\begin{align}
& 2 \pi \int_{-R}^{R} \left( \int_0^{\sqrt{R^2-z^2}} \frac{rdr}{\sqrt{r^2 + (z+p)^2}} \right) dz\\
= & 2 \pi \int_{-R}^{R} \left[\sqrt{r^2+(z+p)^2}\right]_0^{\sqrt{R^2-z^2}} dz\\
= & 2 \pi \int_{-R}^{R} \left(\sqrt{R^2+2zp+p^2} - (z+p)\right) dz\\
= & 2 \pi \left[ \frac{1}{3p} \sqrt{R^2+2zp+p^2}^3 - (\frac{z^2}{2} + p z)\right]_{-R}^R\\
= & 2 \pi \left( \frac{1}{3p}\left( (p+R)^3 - (p-R)^3 \right) - 2pR \right)\\
= & \frac{4\pi R^3}{3\sqrt{a^2+b^2+c^2}}
\end{align}$$
