Calculating in closed form $\int _0^1\int _0^1\frac{1}{1+x y (x+y)} \ dx \ dy$ Integrating with respect to a variable and then to the other one, things look pretty
complicated, but I'm sure you have ideas that might simplify the job to do here.
This time we're talking about 
$$\int _0^1\int _0^1\frac{1}{1+x y (x+y)} \ dx \ dy$$

The bounty moment: after 2 years and 8 months from the releasing moment of the question, it's time for a 300 points bounty for finding the simplest closed-form of the integral!

Supplementary question: Here is an extension of the question for those with a need for more challenging questions.
Calculate
$$\int _0^1\int _0^1\cdots\int _0^1\frac{1}{1+x_1 x_2\cdots x_n (x_1+x_2+\cdots +x_n)} \ \textrm{d}x_1 \ \textrm{d}x_2\cdots \textrm{d}x_n, \ n\ge 2.$$
Last but not least, special greetings to Cleo!

Starting in the polar coordinates,
\begin{align}
&J=\int\limits_0^1\int\limits_0^1 \dfrac{dx\,dy}{1+xy(x+y)} = 2\int\limits_0^{\dfrac\pi4}\int\limits_0^\dfrac1{\cos\varphi}\dfrac{\rho\,d\rho\,d\varphi}{1+\rho^3\cos\varphi\sin\varphi(\cos\varphi+\sin\varphi)} = 2\int\limits_0^{\pi/4}F(\rho,\varphi)\,d\varphi,\\
\end{align}
and substitution
\begin{align}
&\rho=\dfrac1{v \cos\varphi}\tag1
\end{align}
gives
\begin{align}
&F(\rho,\varphi) = \int\limits_1^\infty\dfrac1 {1+\dfrac1{v^3\cos^3\varphi}\cos\varphi\sin\varphi(\cos\varphi+\sin\varphi)}\,\dfrac{dv}{v^3\cos^2\varphi} =\\ &\dfrac1{\cos^2\varphi}\int\limits_1^\infty\dfrac{dv}{\tan\varphi(\tan\varphi+1)+ v^3},\quad\text{or}\\
&J=2\int\limits_1^\infty G(v,\varphi)\, dv,\quad\text{where}\\
&G(v,\varphi) = \int\limits_0^{\pi/4}\dfrac{dv}{\tan\varphi(\tan\varphi+1)+v^3}\cdot\dfrac{d\varphi}{\cos^2\varphi}=\int\limits_0^1\dfrac{dt}{t^2+t+v^3} = \\
&\dfrac1{\sqrt{v^3-\dfrac14}}\left.\arctan\dfrac{t+\dfrac12}{\sqrt{v^3-\dfrac14}}\right|_{t=0}^1 = \dfrac1{\sqrt{v^3-\dfrac14}}\left(\arctan\dfrac{\dfrac32}{\sqrt{v^3-\dfrac14}}-\arctan\dfrac{\dfrac12}{\sqrt{v^3-\dfrac14}}\right),\tag2\\[4pt]
&G(v,\varphi) = \dfrac1{\sqrt{v^3-\dfrac14}}\arctan\dfrac{\dfrac1{\sqrt{v^3-\dfrac14}}}{1+\dfrac3{4\left(v^3-\dfrac14\right)}} = \dfrac2{\sqrt{4v^3-1}}\arctan\dfrac{2\sqrt{4v^3-1}}{4v^3+2}.\tag3
\end{align}
Therefore,
\begin{align}
&J=\int\limits_1^\infty \dfrac4{\sqrt{4v^3-1}}\arctan\dfrac{2\sqrt{4v^3-1}}{4v^3+2}\, dv\approx 0.798965\tag4
\end{align}
(see also Wolfram Alpha).
Substitution
$$u=\sqrt{4v^3-1}$$
leads to
\begin{align}
&v=\dfrac{\sqrt[3]2}2\sqrt[3]{u^2+1},\quad dv=\dfrac{2\sqrt[3]2}3\dfrac1{\sqrt[3]{(u^2+1)^2}}u\,du,\\
&J=\dfrac{4\sqrt[3]2}3\int\limits_\sqrt3^\infty\arctan\dfrac {2u}{3+u^2}\dfrac{du}{\sqrt[3]{(u^2+1)^2}}\tag5\\
&J=\dfrac{4\sqrt[3]2}3\int\limits_\sqrt3^\infty\arctan\dfrac {u-\dfrac13u}{1+u\cdot\dfrac13u}\dfrac{du}{\sqrt[3]{(u^2+1)^2}},\\
&J=\dfrac{4\sqrt[3]2}3\int\limits_\sqrt3^\infty\left(\arctan u - \arctan\dfrac u3\right)\dfrac{du}{\sqrt[3]{(u^2+1)^2}} = J_1-J_2, \tag6\\
\end{align}
(see also Wolfram Alpha), wherein
$$J_1=\dfrac{4\sqrt[3]2}3\int\limits_{\pi/3}^{\pi/2}\dfrac{z\,dz}{\sqrt[3]{\cos^2z}}\\$$
has the closed form.
But the further attempts to get a closed form weren't successful. 
