How to evaluate $S=\int_{\partial D}{\frac{1}{|x-y|^2}}ds(y)$? Let $D=B(0,r_0), x,y \in \mathbb{R}^3$ and $x\in D$.
\begin{align*}
   & S=\int_{\partial D}{\frac{1}{|x-y|^2}}ds(y)=\int_{\partial D}{\frac{1}{(x_1-y_1)^2+(x_2-y_2)^2+(x_3-y_3)^2}}ds(y)
\end{align*}
I want to prove that $S\leq k$ ($k>0$ and $k$ only depends on $r_0$). 
In case $x \ne 0$, my idea is 
\begin{split}
   S^+ &=\int_{y_1^2+y_2^2\leq r_0^2}{\frac{1}{(x_1-y_1)^2+(x_2-y_2)^2+\Big(x_3-\sqrt{r_0^2-y_1^2-y_2^2}\Big)^2}}\frac{r_0}{\sqrt{r_0^2-y_1^2-y_2^2}}dy_1dy_2\\
   &\leq \int_{y_1^2+y_2^2\leq r_0^2}{\frac{1}{\Big(x_3-\sqrt{r_0^2-y_1^2-y_2^2}\Big)^2}}\frac{r_0}{\sqrt{r_0^2-y_1^2-y_2^2}}dy_1dy_2\: \text{ (assume } x_3 \ne 0)
  \end{split}
Let $y_1=r\cos\phi, y_2=r\sin\phi, 0\leq \phi \leq 2\pi$
\begin{split}
   S^+&=\int_0^{2\pi}d\phi\int_0^{r_0}\frac{1}{\Big(x_3-\sqrt{r_0^2-r^2}\Big)^2}\frac{r_0}{\sqrt{r_0^2-r^2}}rdr\\
&=2r_0\pi\left(\frac{1}{-x_3}-\frac{1}{r_0-x_3}\right)\\
  \end{split}
with $-r_0<x_3<r_0, x_3 \ne 0$, I can't show that $\left(\frac{1}{-x_3}-\frac{1}{r_0-x_3}\right)$ is bounded. Please help me, thank you so much!
 A: You may use spherical coordinates and the cosine-theorem:
$$
\frac{1}{|x-y|^2} = \frac{1}{|x|^2+|y|^2-2|x||y|\cos(\theta_{x,y})}
$$
where $\theta_{x,y}$ is the angle between $x$ and $y$.
Putting it into the integral you have:
$$
\int\limits_{|y|= r_0}\frac{ds}{|x-y|^2} = 
\int\limits_{|y|= r_0}\frac{ r_0^2\sin(\theta) d\phi d\theta}{|x|^2+r_0^2-2|x|r_0\cos(\theta)}=2\pi r_0^2\int\limits_{-1}^1\frac{ d\cos(\theta)}{|x|^2+r_0^2-2|x|r_0\cos(\theta)}
$$
where we chose to allign $x$ along the z-axis, which we are free to do due to symmetry. Then $\theta_{x,y}$ simply becomes the coordinate $\theta$.
A: Clearly, 
$$
F(y)=\int_{|y|=r_0}\frac{ds}{|s-y|^2}
$$
is a function of $|y|$. For simplicity, assume that $y=(|y|,0,0)$.
Since the sphere is a solid by revolution, and in particular of revolution of the graph of the function $f(x)=\sqrt{r^2-x^2}$ around the $x-$axis, then 
$$
F(y)=2\pi\int_{-r_0}^{r_0} f(x)\sqrt{1+\big(f'(x)\big)^2}g(x)\,dx 
=\cdots=2\pi r_0\int_{-r_0}^{r_0} g(x)\,dx
$$
where 
$$
g(x)=\frac{1}{\big|\big(x,\sqrt{r_0^2-x^2},0\big)-(|y|,0,0)\big|^2}=\frac{1}{(x-|y|)^2+\big(\sqrt{r_0^2-x^2}\big)^2}=\frac{1}{r_0^2-2x|y|+|y|^2}
$$
and hence
$$
F(y)=2\pi r_0\int_{-r_0}^{r_0}\frac{dx}{r_0^2+|y|^2-2|y|x}=
\left.-\frac{\log(r_0^2+|y|^2-2|y|x)}{2|y|}\right|_{-r_0}^{r_0}
=\frac{\log\left(\frac{r^2_0+2r_0|y|+|y|^2}{r^2_0-2r_0|y|+|y|^2}\right)}{2|y|}=
\frac{
\log\left(
\frac{r_0+|y|}{r_0-|y|}
\right)
}{|y|}
$$
Clearly,
$$
\lim_{|y|\to r_0^-}F(y)=\infty
$$
