What substitution would I make to integrate this? Problem:

The integral is gonna be:
$$\int^{1}_{0}\int^{1}_{0}4xy\sqrt{x^2+y^2} dy \, dx$$
But I'm quite rusty with my calculus, and this is mainly for a statistics course. I know that $x$ and $y$ are bounded within a unit square region in the $xy$-plane, and I was considering making some sort of $\cos$ or $\sin$ substitution to get rid of that pesky square root, but then I'm not sure what limits to use, and I suspect the $dy dx$ would just become $d\theta d\theta$, which doesn't make sense.
Any help or guidance is appreciated. Please try to keep your explanation simple, and don't assume I understand things. It's been a while since my last calculus course :) Thank you!
 A: \begin{align}
  E(Z) &= \int_{0}^{1} \int_{0}^{1} Z f(x,y) \, dx \, dy \\
  &= \int_{0}^{1} \int_{0}^{1}  4xy\sqrt{x^2+y^2} \, dx \, dy \\
  &= \int_{0}^{1} 2y \left( \int_{0}^{1}  2x\sqrt{x^2+y^2} \, dx \right) dy \\
  &= \int_{0}^{1} 2y
     \left(
       \int_{0}^{1}  \sqrt{u+y^2} \, du
     \right)
      dy \tag{$u=x^2$} \\
  &= \int_{0}^{1} 2y
     \left[
       \frac{2}{3} (u+y^{2})^{3/2}
     \right]_{u=0}^{1} dy \\
  &= \int_{0}^{1} \frac{4y}{3}
     \left[
       (1+y^{2})^{3/2}-y^{3}
     \right] dy \\
  &= \frac{4}{15}
     \left[
       (1+y^{2})^{5/2}-y^{5}
     \right]_{y=0}^{1} \\
  &= \frac{8}{15}(2\sqrt{2}-1)
\end{align}
A: $$\int^{1}_{0}\int^{1}_{0}4xy\sqrt{x^2+y^2} dy \, dx=$$
$$\int^{\pi/4}_{0}\int^{\sec\theta}_{0}4r^4\sin\theta \cos \theta drd\theta+\int^{\pi/2}_{\pi/4}\int^{cosec \theta}_{0}4r^4\sin\theta \cos \theta drd\theta $$
Explanation:-
$$x=r\cos \theta$$
$$y=r\sin \theta$$
$$dxdy=rdrd\theta$$

$0\leq\theta \leq \pi/4$ and for the variation in $r$,$r:0\to sec \theta(\because x=1$  can be translated as$ r\cos\theta=1)$. Simillarly other region.
A: An alternative approach:
$$ \begin{eqnarray*}\iint_{(0,1)^2}4xy\sqrt{x^2+y^2}\,dx\,dy&=&\iint_{(0,1)^2}\sqrt{X+Y}\,dX\,dY\\&=&2\iint_{0\leq Y\leq X\leq 1}\sqrt{X+Y}\,dX\,dY\\&=&2\iint_{(0,1)^2}X\sqrt{X}\sqrt{K+1}\,dX\,dK\\&=&2\int_{0}^{1}X\sqrt{X}\,dX\int_{0}^{1}\sqrt{K+1}\,dK\end{eqnarray*}$$
by exploting the substitution $x=\sqrt{X}, y=\sqrt{Y}$, then symmetry, then the substitution $Y=KX$, then Fubini's theorem. The last integrals are elementary and they respectively equal $\frac{2}{5}$ and $\frac{2}{3}(2\sqrt{2}-1)$.
