An integral related to the distance distribution between points inside the unit disk At the moment I am interested in some applications of geometric probability theory, like finding the average distance $\langle D\rangle$ of two points randomly picked from the unit disk. It is not too difficult to set up an integral which correctly describes this distance:
Pick two points $X_1,X_2$ inside the unit disk (one of them can, by symmetry, chosen to be located on the y-axis) and denote by $r_j=|X_j|$ their distance from the origin. 
Then $r_j$ are indenpendently distributed random variables with probability distribution $dp_j=2r_j$. Furthermore, the angle between these two points $\phi$ is uniformly distributed with pd $dp_{\phi}=1/\pi$.
From elementary geometry (law of cosines) is easy to see that $D^2=r_1^2+r_2^2-2 r_1r_2\cos(\phi)$ and therefore ($A$ is a cylinder with $r=1/2, z=1$   )
$$
\langle D\rangle=\int_ADdp_1dp_2dp_{\phi}=
\\\frac{4}{\pi}\int_0^1\int_0^1\int_0^{\pi}r_1 r_2\sqrt{r_1^2+r_2^2-2 r_1r_2\cos(\phi)}dr_1dr_2d\phi
$$
After struggeling for a few hours, i was googeling around to find out that 

$$\langle D\rangle=\frac{128}{45\pi} \quad\quad (\star)$$

Fighting a bit longer, i am now ready to surrender and asking you 

How on Earth/Mars/Beteigeuze $(\star)$ can be proven without plodding through a swamp of nightmarish antiderivatives?

I tried a lot, from bruteforce methods over various changes of coordinate systems to series expansions but everything ended up in a big mess (interestingly for 3 dimensions the calculations go through smoothly)...
Edit:
Just in case somebody needs a bit of visual input, i created a Desmos chart

Edit2:
 A: Due to a powerful theorem due to Crofton it is possible to proof $(\star)$  by indirect means. The theorem states for our  particular problem (we introduce variable disk radius $R$)
$\frac{d\langle D\rangle}{dR}=2(\langle D\rangle_1-\langle D\rangle)\frac{2}{R} \quad (\star\star)$ 
Here $\langle D\rangle_1$ is the average distance if one of the points sits at the circumference of the disk
($2/R$ comes from the quotient of the infinitesimal change of volume and the disk volume)
which is easy enough to calculate since this distance is just the average length of the inscribed arc length of a second circle intersecting our disk at some distance $y$ away from our point at the circumference and is given by $2y\Theta(y)$ where $\Theta(y)=\arccos(y/2R)$ is the angle betweeen the Diameter and the line connecting the two points on the circle
$$
\langle D\rangle_1= \frac{2}{\pi R^2} \int_0^{2R}yy\arccos(y/2R)dy=\frac{32}{9\pi}
$$
Putting this back into $(\star\star)$ and solving the ODE demanding $\langle D\rangle_{R=0}=0$ yields

$$
\langle D\rangle=\frac{128}{45\pi R}
$$

QED
The direct Evaluation of the integral in my Question nevertheless stays a mystery
Some visualisation:

