I want to evaluate$$\int_{S^2} \frac{dS}{((x-a)^2 +y^2+z^2)^{1/2}}$$ where $a >1$ and $S$ is the unit sphere.
I'm not sure how to do this using only multivariable calculus techniques. My only idea was to use the fact that the function $(x^2+y^2+z^2)^{-1/2}$ is harmonic away from the origin and use the mean value formula.
We can use a few symmetry conditions to make our lives easier. Note that your proposed integral is equivalent to
$$\iint_{(x+a)^2+y^2+z^2=1} \frac{1}{\sqrt{x^2+y^2+z^2}}\:dS = \iint_{(x+a)^2+y^2+z^2=1} \frac{(x+a,y,z)\cdot(x+a,y,z)}{\sqrt{x^2+y^2+z^2}}\:dS$$
Thus we can use the divergence theorem
$$= \iint_{(x+a)^2+y^2+z^2\leq 1} \frac{2}{\sqrt{x^2+y^2+z^2}}-\frac{ax}{(x^2+y^2+z^2)^{\frac{3}{2}}}\:dV$$
In spherical coordinates centered around the $x$ axis, the boundary can be rewritten as
$$x^2+y^2+z^2+2ax+a^2-1 = 0 \implies \rho^2+2a\rho\cos\phi + a^2-1 = 0$$
The best way to set up the integral is with $\phi$ first like so:
$$\int_0^{2\pi} \int_{a-1}^{a+1} \int_{\cos^{-1}\left(\frac{1-a^2}{2a\rho}-\frac{\rho}{2a}\right)}^\pi 2\rho\sin\phi - a\cos\phi\sin\phi \:d\phi\:d\rho\:d\theta$$
$$ = 2\pi\int_{a-1}^{a+1} 2\rho + \frac{5-3a^2}{4a} - \frac{9\rho^2}{8a} -\frac{(1-a^2)^2}{8a\rho^2}\:d\rho$$
$$= 2\pi \left[4a + \frac{5-3a^2}{2a}-\frac{3}{8a}(6a^2+2) - \frac{a^2-1}{4a}\right] = \frac{4\pi}{a}$$
where you can see the answer was more easily retrieved by the Mean Value Property of harmonic functions (i.e. $\int_{S^2} f = 4\pi f(0)$ due to the symmetries in the problem), but this is a purely multivariable calculus integration way of proving the result without resorting to PDEs.