Is there a nice way to evaluate $\iiint_{E}\, \frac{dx\,dy \,dz}{\sqrt{x^2+y^2+(z-b)^2}}$ where $E:x^2+y^2+z^2\leq a^2$ 
Is there a nice way to evaluate $$\displaystyle\iiint_{E}\, \dfrac{dx\,dy\,dz}{\sqrt{x^2+y^2+(z-b)^2}}$$ where   $E:x^2+y^2+z^2\leq a^2$ with
    $0<a<b$

If I use the standard spherical coordinates (Is there a better transformation?)
$$x=p\sin(\phi)\cos(\theta) $$
$$y=p\sin(\phi)\sin(\theta)$$
$$z=p\cos(\phi)$$
$$|J|=p^2\sin(\phi)$$
$$\iiint_{E}\, \frac{dx\,dy \,dz}{\sqrt{x^2+y^2+(z-b)^2}} $$
$$=\int_{\phi=0}^{\pi}\,\int_{\theta=0}^{2\pi}\,\int_{p=0}^{a}\, \frac{p^2\sin(\phi)\,dp\,d\theta \,d\phi}{\sqrt{p^2\sin^2(\phi) + (p\cos(\phi)-b)^2}}$$ 
$$=\int_{\phi=0}^{\pi}\,\int_{\theta=0}^{2\pi}\,\int_{p=0}^{a}\, \frac{p^2\sin(\phi)\,dp\,d\theta \,d\phi}{\sqrt{(p-b\cos(\phi))^2 + b^2\sin^2(\phi)}}$$
Which I'm finding bit difficult. 
 A: Fubini split w.r.t. the $z$-component. Let $J$ be the integral to be calculated. I will use 
$B(0,r)$ for a ball centered in $0$ with radius $r$ (in the dimension that 
can be extracted from the context.) 
Then:
$$
\begin{aligned}
J &=
\iiint_{B(0,a)}\frac 1{\sqrt{x^2+y^2+(b-z)^2}}\; dx\; dy\; dz
\\
&=
\int_{z\in B(0,a)}dz
\iint_{(x,y)\in B(0,\sqrt{a^2-z^2})}\frac 1{\sqrt{x^2+y^2+(b-z)^2}}\; dx\; dy
\\
&=
\int_{-a}^a dz
\int_{r\in[0,\sqrt{a^2-z^2}]}
\int_{t\in[0,2\pi]}
\frac 1{\sqrt{r^2+(b-z)^2}}\; r\; dr\; dt
\\
&=
\int_{-a}^a dz
\int_0^{\sqrt{a^2-z^2}}
2\pi\Big(\ (r^2+(b-z)^2)^{1/2}\ \Big)'\; dr
\\
&=
2\pi
\int_{-a}^a dz\;
\Big[\ (r^2+(b-z)^2)^{1/2}\ \Big]_0^{\sqrt{a^2-z^2}}
\\
&=
2\pi
\int_{-a}^a dz\;
\Big[\ (a^2-z^2+(b-z)^2)^{1/2}-(0^2+(b-z)^2)^{1/2}\ \Big]
\\
&=
2\pi
\int_{-a}^a dz\;
\Big[\ (a^2-2bz+b^2)^{1/2}-(b-z)\ \Big]
\\
&=
2\pi
\left[\ -\frac 1{2b}\cdot \frac 1{3/2}(a^2-2bz+b^2)^{3/2}
+\frac 12(b-z)^2\ \right]_{-a}^a
\\
&=
\frac {2\pi}{3b}\Big[\ (b+a)^3-(b-a)^3\ \Big]
-
\pi\Big[\ (b+a)^2-(b-a)^2\ \Big]
\\
&
=
\frac {2\pi}{3b}\cdot 2a^3
=
\frac {4\pi\, a^3}{3b}\ .
\end{aligned}
$$
Check using sage:
sage: var( 'r,t,z,a,b' );
sage: assume(a>0)
sage: assume(b>a)
sage: J = integral( integral( integral( r/sqrt(r^2+(z-b)^2), t,0,2*pi ), r,0,sqrt(a^2-z^2) ), z,-a,a )
sage: J.canonicalize_radical()
4/3*pi*a^3/b


Here is the computation of the integral for its version "without the square root in the denominator"...
$$
\begin{aligned}
J &=
\iiint_{B(0,a)}\frac 1{x^2+y^2+(b-z)^2}\; dx\; dy\; dz
\\
&=
\int_{z\in B(0,a)}dz
\iint_{(x,y)\in B(0,\sqrt{a^2-z^2})}\frac 1{x^2+y^2+(b-z)^2}\; dx\; dy
\\
&=
\int_{-a}^a dz
\int_{r\in[0,\sqrt{a^2-z^2}]}
\int_{t\in[0,2\pi]}
\frac 1{r^2+(b-z)^2}\; r\; dr\; dt
\\
&=
\int_{-a}^a dz
\int_0^{\sqrt{a^2-z^2}}
2\pi\cdot\frac 12\Big(\ \ln(r^2+(b-z)^2)\ \Big)'\; dr
\\
&=
\pi
\int_{-a}^a dz\;
\Big[\ \ln(r^2+(b-z)^2)\ \Big]_0^{\sqrt{a^2-z^2}}
\\
&=
\pi
\int_{-a}^a dz\;
\Big[\ \ln(a^2-z^2+(b-z)^2)-\ln(0^2+(b-z)^2)\ \Big]
\\
&=
\pi
\int_{-a}^a dz\;
\Big[\ \ln(a^2-2bz+b^2)-2\ln(b-z)\ \Big]
\\
&=
\pi
\left[\ -\frac 1{2b}(a^2-2bz+b^2)(\ln(a^2-2bz+b^2)-1)+2(b-z)(\ln(b-z)-1)\ \right]
_{-a}^a
\\
&=
\pi
\left[\ 
-\frac 1{2b}(b-a)^2(2\ln(b-a)-1)
+\frac 1{2b}(b+a)^2(2\ln(b+a)-1)
\right.
\\
&\qquad\qquad
\left.
+2(b-a)(\ln(b-a)-1)
-2(b+a)(\ln(b+a)-1)
\ \right]
\\
&=
2\pi\ln(b+a)\cdot\left[\frac 1b(b+a)^2-2(b+a)\right]
-
2\pi\ln(b-a)\cdot\left[\frac 1b(b-a)^2-2(b-a)\right]
+2\pi a
\\
&=
2\pi\frac 1b(b^2-a^2)\ln\frac {b-a}{b+a}+2a\pi\ .
\end{aligned}
$$
A: I would try cylindrical coordinates:
$$
\begin{cases}
x=\rho\cos\phi\\
y=\rho\sin\phi\\
z=z
\end{cases}
\implies
\begin{cases}
x^2+y^2=\rho^2\\
dx\,dy=\rho\,d\rho\,d\phi
\end{cases}
$$
Which should give you:
\begin{align}
\iiint_E \frac{\rho}{\rho^2+(z-b)^2}d\phi\,d\rho\,dz&=\iint_E \frac{2\pi\rho}{\rho^2+(z-b)^2}d\rho\,dz\\
&=\int_{-a}^a\int_0^{\sqrt{a^2-z^2}} \frac{2\pi\rho}{\rho^2+(z-b)^2}d\rho\,dz\\
\end{align}
To me, it seems a little easier now.
A: Starting with your final integral, do four things:
1) Since the integrand does not depend on $\theta$, perform the $\theta$ integration, which puts a factor of $2\pi$ in front.
2) Multiply out and simplify the things in the denominator.
3) Make the change of variables $p = b\, r$.
4) Make the change of variables $u = \cos\phi$.
After all that, we are left with:
$$
2\pi b^2 \int_0^{a/b}dr\, r^2 \int_{-1}^1 du\, \frac{1}{\sqrt{1 - 2 r u + r^2}}
$$
Note that the quantity $1/\sqrt{1 - 2 r u + r^2}$ is the generating function for Legendre Polynomials:
$$
\frac{1}{\sqrt{1 - 2 r u + r^2}}
\;=\; \sum_{n=0}^{\infty} P_n(u)\, r^n
$$
Also note that
\begin{equation}
\int_{-1}^1 du\, P_n(u) \;=\;
\begin{cases}
2 & n = 0\\
0 & n\ne 0
\end{cases}
\end{equation}
Putting this all together:
\begin{align}
2\pi b^2 \int_0^{a/b}dr\, r^2 \int_{-1}^1 du\, \frac{1}{\sqrt{1 - 2 r u + r^2}}
&=
2\pi b^2 \int_0^{a/b}dr\, r^2 \int_{-1}^1 du\,\sum_{n=0}^{\infty} P_n(u)\, r^n\\
&=
4\pi b^2 \int_0^{a/b}dr\, r^2\\
&= \frac{4\pi}{3}\,\frac{a^3}{b}
\end{align}
I like to check things numerically in case of mistakes, and having done so in Mathematica I find that this expression matches the numerical integration of your final integral for all the values of $b \ge a$ that I tried.
A: An easier way to slice things up is to cut the sphere $x^2+y^2+z^2\leq a^2$ into spherical shells centered at $(0,0,b)$.
First of all, observe that the desired result is the integral over the sphere $x^2+y^2+z^2\leq a^2$ of the reciprocal distance to the point $(0,0,b)$. This is the motivation for slicing things up into shells centered at $(0,0,b)$.
In the picture below the sphere $x^2+y^2+z^2\leq a^2$ is pink, and a few of the shells over which we’ll integrate are indicated by the grey arcs (and one heavier black arc) within the pink area.

The point of this approach is that the reciprocal distance to $(0,0,b)$ is constant on each shell.
Parameterize the shells by the $z$-coordinate at which they cross the blue line below. The $z$-values run from $-a$ to $a$. The points of the “shell at $z$” are at distance $b-z$ from $(0,0,b)$, so the answer to the original question is
$$\int_{-a}^{a}{\frac1{b-z}}\left(\mbox{ surface area of shell at $z$ }\right)\,dz$$
The shell at $z$ is a “spherical cap,” and it’s been known since Archimedes’ time that the surface area of a spherical cap of height $h$ cut from a sphere of radius $r$ is $2\pi r h$. A bit of geometry and algebra shows that the spherical cap that is this problem’s “shell at $z$” has height $\dfrac{a^2-z^2}{2b}$, and as previously observed, is part of a sphere of radius $b-z$. Therefore the answer to the original question is
$$\int_{-a}^{a}{\frac1{b-z}}2\pi(b-z)\left(\dfrac{a^2-z^2}{2b}\right)\,dz=\frac{\pi}b\int_{-a}^{a}\left({a^2-z^2}\right)\,dz=\frac{4\pi a^3}{3b}$$
