I saw a meme about $\int_{0}^1 \int_0^1\int_0^1 \sqrt{x^2 + y^2 + z^2} \space dx dy dz$ on Facebook (Wolfram Alpha only giving an approximate answer) and the closed form of the answer is supposed to be $\frac{\pi}{24} + \frac{1}{4} \left( \sqrt{3} + \log(7 + 4\sqrt{3}) \right)$, which I verified with Mathematica, but Mathematica isn't offering any steps or hints as to how it arrived at that answer.
I've tried it a few times myself and in cylindrical coordinates end up with $2\int_0^1 \int_0^{\frac{\pi}{4}}\int_0^{\sec\theta} r\sqrt{r^2 + z^2} dr d\theta dz = \frac{2}{3}\int_0^1\int_0^{\frac{\pi}{4}}\left( \sec^2\theta + z^2\right)^{\frac{3}{2}} - z^3 d\theta dz$ by just considering the first half of the cube swept out and doubling it. The $d\theta$ integral was looking pretty awful from what Wolfram had to offer. I then tried to do a $z = \sec\theta \tan u $ substitution (treating $\sec \theta $ as a constant) after swapping the order of integration. This leads to having to evaluate $\int_0^{\arctan\left( cos \theta\right)}\sec^5 u \space du$ and that itself involves evaluating a number of terms such as $\sin 3u$ at these limits.
I hope I've made a simple mistake along the line but can anyone spot a nicer/simpler approach?
Cheers.