Evaluate $\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{(z-w)^2+(x-y)^2} \, dw \, dz \, dy \, dx$

From this page I found an interesting equality:$$\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{(z-w)^2+(x-y)^2} \, dw \, dz \, dy \, dx=\frac{1}{15} \left(\sqrt{2}+2+5 \log \left(\sqrt{2}+1\right)\right)$$ It's known as the average distance between two random points in the unit square, and higher dimensional results are also obtained in that page. My question is: is there any elementary method to obtain this result? Thank you. $$$$Update: In comments one user offered a link which solves the problem by basic methods. Still, alternative methods for this problem (as well as 3 dimension's analogue, see the link above) are welcomed.

from Hot Weekly Questions - Mathematics Stack Exchange