Aquí es una manera alternativa que evita su integral.
Escrito el integrando en términos de un geométrica de la suma tenemos:
$$\frac{1}{1 + y(x^2 - x)} = \sum_{n = 0}^\infty y^n (x - x^2)^n, \quad |x|, |y| < 1,$$
entonces
\begin{align}
I &= \int_0^1 \int_0^1 \frac{1}{1 + y(x^2 - x)} \, dy dx\\
&= \sum_{n = 0}^\infty \int_0^1 y^n \, dy \int_0^1 x^n (1 - x)^n \, dx\\
&= \sum_{n = 0}^\infty \frac{1}{n + 1} \cdot \operatorname{B}(n + 1, n + 1)\\
&= \sum_{n = 0}^\infty \frac{1}{n + 1} \cdot \frac{(n!)^2}{(2n + 1)!}\\
&= \sum_{n = 0}^\infty \frac{1}{(n + 1)(2n + 1) \binom{2n}{n}}.\tag1
\end{align}
Aquí $\operatorname{B}(x,x)$ denota la función Beta, mientras que $\binom{2n}{n} = \frac{(2n)!}{(n!)^2}$ es la central de coeficiente binomial.
Para encontrar la suma en (1), el siguiente (razonablemente?) conocida la serie de Maclaurin se pueden recuperar (para una prueba, ver aquí)
$$\frac{\sin^{-1} x}{\sqrt{1 - x^2}} = \sum_{n = 0}^\infty \frac{2^{2n} x^{2n + 1}}{(2n + 1) \binom{2n}{n}}, \quad |x| < 1.$$
El cumplimiento de una sustitución de $x \mapsto \sqrt{x}$, después de la reorganización de uno ha
$$\frac{\sin^{-1} \sqrt{x}}{\sqrt{x} \sqrt{1 - x}} = \sum_{n = 0}^\infty \frac{2^{2n} x^n}{(2n + 1) \binom{2n}{n}}, \quad |x| < 1.$$
La integración de ambos lados con respecto a $x$ de $0$ a $\frac{1}{4}$ conduce a
$$\sum_{n = 0}^\infty \frac{1}{(n + 1)(2n + 1) \binom{2n}{n}} = 4 \int_0^{\frac{1}{4}} \frac{\sin^{-1} \sqrt{x}}{\sqrt{x} \sqrt{1 - x}} \, dx.$$
Una sustitución de $x = \sin^2 u$ en la integral fácilmente conduce a
\begin{align}
\sum_{n = 0}^\infty \frac{1}{(n + 1)(2n + 1) \binom{2n}{n}} &= 8 \int^{\frac{\pi}{6}}_0 u \, du = \frac{\pi^2}{9},
\end{align}
a partir de la cual se deduce que:
$$\int_0^1 \int_0^1 \frac{1}{1 + y(x^2 - x)} \, dy dx = \frac{\pi^2}{9}.$$