I know that the average value of a function $f(x,y,z)$ over a volume $D$ is given by
$$ \bar{f} = \frac{1}{V} \int\int\int_D f(x,y,z) dV $$
I'm trying to set up the integral of the average distance from a point some distance a from the origin of a hollow sphere to the surface of the sphere (image from here):

I used this equation to find $d$: $R^2=a^2+d^2+2adcos(\theta)$ $$ \Rightarrow d = \sqrt{R^2-a^2sin^2(\theta)}-acos(\theta) $$ My problem is I am now having difficulty setting up the integral. I know I'm supposed to get: $$ \bar{r} = \frac{3R}{2} \int_0^1 \int_0^{\pi} \sqrt{1-k^2sin^2(\theta)}k^2sin(\theta)dkd\theta $$ About the closest I can get to that is: $$ \bar{d} = \frac{3}{4\pi R^3} \int_0^1 \int_0^{\pi} \int_0^{2\pi} \bigg(\sqrt{R^2-a^2sin^2(\theta)}-acos(\theta)\bigg)k^2sin(\theta)dkd\phi d\theta $$ I'm clearly conceptually missing something, so any advice for what I'm misunderstanding or hints on how to set up this integral would be greatly appreciated.
If one want to average over the hollow sphere instead of average over all possible directions with respect to the reference point, one should use the green $\theta$ instead of blue $\theta$. They are two different angles and expressions will become very complicated if one use the wrong angle for the averaging scheme it is not intended for.
To avoid confusion with the differential operator, I will use $r$ instead of $d$ to denote the distance between the reference point and point on the hollow sphere. Based on the given picture, I will also assume the reference point is inside the hollow sphere.
Case I - average over the hollow sphere.
In this case, the $\theta$ is the green $\theta$. We have
$$r = \sqrt{a^2 + R^2 - 2aR\cos\theta}$$
The integral I will evaluate is
$$\bar{r} \stackrel{def}{=} \frac{1}{4\pi R^2}\int_{S^2} r R^2 d\Omega = \frac{1}{4\pi}\int_{S^2} r d\Omega = \frac{1}{4\pi}\int_{S^2} \frac{a^2+R^2-2aR\cos\theta}{r} d\Omega $$ where $d\Omega = \sin\theta d\theta d\phi$ is the area element of unit sphere $S^2$.
To compute the integral, we use something called Multipole expansion. In general, if we have two distinct points $\vec{p}$ and $\vec{q}$ in $\mathbb{R}^3$, we have following expansion: $$\frac{1}{|\vec{p} - \vec{q}|} = \sum_{\ell = 0}^\infty \frac{r_<^\ell}{r_>^{\ell+1}} P_\ell(\cos\alpha) \quad\text{ where }\quad \begin{cases} r_< &= \min\{ |\vec{p}|, |\vec{q}| \},\\ r_> &= \max\{ |\vec{p}|, |\vec{q}| \},\\ \alpha &= \cos^{-1}\frac{\vec{p}\cdot\vec{q}}{|\vec{p}||\vec{q}|} \end{cases} $$ and $P_\ell(x)$ is the Legendre polynomial of order $\ell$.
The key properties we need is $P_\ell(\cos\theta)$ forms a orthogonal family of polynomials over $S^2$.
More precisely,
$$\frac{1}{4\pi} \int_{S^2} P_\ell(\cos\theta)P_{\ell'}(\cos\theta) d\Omega = \begin{cases} \frac{1}{2\ell+1},& \ell = \ell',\\ 0, & \ell \ne \ell' \end{cases}$$ Notice $P_0(x) = 1$, $P_1(x) = x$, we have
$$\begin{align}\bar{r} &= \frac{1}{4\pi}\int_{S^2} \left((a^2+R^2)P_0(\cos\theta) - 2aRP_1(\cos\theta)\right)\left( \sum_{\ell=0}^\infty \frac{a^\ell}{R^{\ell+1}}P_\ell(\cos\theta) \right) d\Omega\\ &= (a^2+R^2)\left(\frac{1}{R}\right) - (2aR)\left(\frac{a}{3R^2}\right)\\ &= R + \frac{a^2}{3R} \end{align} $$
Case II - all directions with respect to reference point is equally likely.
In this case, the $\theta$ is the blue $\theta$. We have
$$(r\cos\theta + a)^2 + (r\sin\theta)^2 = R^2 \iff (a\cos\theta+r)^2 + (a\sin\theta)^2 = R^2\\ \implies r = \sqrt{R^2 - a^2\sin\theta^2} - a\cos\theta$$ The average we want becomes
$$\bar{r} = \frac{1}{4\pi}\int_{S^2} \left(\sqrt{R^2 - a^2\sin\theta^2} - a\cos\theta \right) d\Omega$$ where $d\Omega = \sin\theta d\theta d\phi$ again but for a different $\theta$.
Notice under the mapping $\theta \mapsto \pi - \theta$, $\cos\theta$ changes sign while the area element $d\Omega$ and the domain of integration remains unchange. the $a \cos\theta$ term doesn't contribute anything. As a result,
$$\bar{r} = \frac{1}{4\pi}\int_0^\pi \int_0^{2\pi} \sqrt{R^2 - a^2\sin\theta^2} \sin\theta d\phi d\theta = \frac12\int_0^\pi \sqrt{R^2 - a^2\sin\theta^2} \sin\theta d\theta $$ Change variable to $t = a\cos\theta$ and notice the intergrand is even under the map $t \mapsto -t$, we have $$\bar{r} = \frac{1}{a}\int_0^a \sqrt{R^2-a^2 + t^2} dt = \frac{R}{2} + \frac{R^2-a^2}{2a}\sinh^{-1}\left(\frac{a}{\sqrt{R^2-a^2}}\right) $$