By symmetry:
$$I=4\int _0^1\int _0^1\int _0^x\int _0^z\sqrt{(z-w)^2+(x-y)^2} \, dw \, dy \, dx dz$$
Substituting $w \to z w, \quad y \to x y$:
$$I=4\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{z^2(1-w)^2+x^2(1-y)^2}~ x z \, dw \, dy \, dx dz$$
Substituting $w \to 1-w, \quad y \to 1-y$:
$$I=4\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{z^2 w^2+x^2 y^2}~ x z \, dw \, dy \, dx dz$$
Substituting $x^2 = u, \quad z^2 = v$:
$$I=\int _0^1\int _0^1\int _0^1\int _0^1\sqrt{v w^2+u y^2}~ \, dw \, dy \, du dv$$
Substituting $v w^2=p, \quad u y^2=q$:
$$I=\int _0^1\int _0^1\int _0^{y^2}\int _0^{w^2} \sqrt{p+q}~ \, dp dq \frac{dy dw}{y^2 w^2}$$
$$I=\frac{2}{3} \int _0^1\int _0^1\int _0^{y^2}\left((q+w^2)^{3/2}-q^{3/2} \right) dq \frac{dy dw}{y^2 w^2}$$
$$I=\frac{4}{15} \int _0^1\int _0^1\left((y^2+w^2)^{5/2}-y^5-w^5 \right) \frac{dy dw}{y^2 w^2}$$
By symmetry:
$$I=\frac{8}{15} \int _0^1\int _0^w\left((y^2+w^2)^{5/2}-y^5-w^5 \right) \frac{dy dw}{y^2 w^2}$$
Substituting: $y = w s$:
$$I=\frac{8}{15} \int _0^1\int _0^1 w^2 \left((1+s^2)^{5/2}-s^5-1 \right) \frac{ds dw}{s^2 }$$
$$I=\frac{8}{45} \int _0^1 \left((1+s^2)^{5/2}-s^5-1 \right) \frac{ds }{s^2 }$$
The last single integral can be evaluated by different methods, and the result agrees with the one from the OP.
One good method here is integration by parts. I leave this as an exercise.
One can exploit a symmetry and split the integral in half along the line $y=x$ then use the following modified polar coordinates:
$$x = s^{\frac{2}{3}}\cos\theta \hspace{10 pt} y = s^{\frac{2}{3}}\sin\theta$$
$$\implies 2\cdot \frac{2}{3}\int_0^{\frac{\pi}{4}} \int_0^{\sec^{\frac{3}{2}}\theta} \frac{\log\left(s^{\frac{4}{3}}\right)}{\sqrt{\cos\theta+\sin\theta}}dsd\theta = \frac{16}{9}\int_0^{\frac{\pi}{4}} \frac{\sec^{\frac{3}{2}}\theta}{\sqrt{\cos\theta+\sin\theta}}\left(\frac{3}{2}\log(\sec\theta)-1\right)d\theta$$
then let $x=\tan\theta$
$$ \implies \frac{4}{3}\int_0^1 \frac{\log(1+x^2)}{\sqrt{1+x}}dx - \frac{16}{9}\int_0^1 \frac{1}{\sqrt{1+x}}dx$$
The integral on the right evaluates to $\frac{32}{9}(\sqrt{2}-1) \equiv \frac{32}{9}a$. The integral on the left becomes
$$ = \frac{8}{3}\sqrt{2}\log 2 - \frac{16}{3}\int_0^1 \frac{x\sqrt{1+x}}{1+x^2}dx$$
We have one and a half of the three terms the user posted, but this last integral is tricky and isn't yielding to a variety of methods. I will try to finish later, but in the meantime if anyone has any clever suggestions for this last integral I will be happy try it.
Best Answer
This is easily tackled by a series of substitution, let $I$ be your integral, then $$\begin{aligned}I &= \int_0^{2\pi} \int_0^{2\pi} \log(3-2\cos \frac y2 \cos(x+\frac y2)-\cos y) dxdy\\ &= \int_0^{2\pi} \int_0^{2\pi} \log(3-2\cos \frac y2 \cos x-\cos y) dxdy\\ &= 4\pi\int_0^{\pi} \log\left[ \frac{1}{2} \left(-\cos \frac{y}{2}+2 \sec \frac{y}{2}+\sqrt{\cos ^2\frac{y}{2}+4 \sec ^2\frac{y}{2}-5}\right)\right] dy\\ &= 8\pi\int_0^1 \log\left(\frac{2-u^2+\sqrt{u^4-5u^2+4}}{2u}\right) \frac{1}{\sqrt{1-u^2}} du \\ &= 64\pi \int_{\pi/6}^{\pi/2} \log(2\sin t) \frac{2-\cos 2 t}{5-4 \cos 2 t} dt \end{aligned}$$ where $4\sin^2 t = 2-u^2+\sqrt{u^4-5u^2+4}$, and $\int_0^1 \frac{\log(2u)}{\sqrt{1-u^2}}du = 0$ is used.
The last integral is straightforward by writing $$\log(2\sin t) = -\sum_{k=1}^\infty \frac{\cos 2kt}{k}\qquad \frac{2-\cos 2 t}{5-4 \cos 2 t} = \frac{1}{2}+\frac{1}{2}\sum_{k=1}^\infty\frac{\cos 2kt}{2^k}$$
Therefore $$I = \sum_{r\geq 0, s>0} \frac{a_{r,s}}{2^r s^2}$$ where $a_{r,s}$ only depends on $r+s \pmod{3}$, the expression of $I$ in terms of polygamma follows immidiately.