You can do what Nate Eldredge suggested. Otherwise, you can use the moment generating function $M$ of the bivariate normal distribution $N(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2,\rho)$ given by 
$$M(t_1,t_2)=\exp(\boldsymbol\mu'\boldsymbol t+\tfrac12\,\boldsymbol t'\Sigma\boldsymbol t)
$$
for $\boldsymbol t=(t_1,t_2)$, $\boldsymbol\mu=(\mu_1,\mu_2)$, 
$\Sigma=\begin{pmatrix}\sigma_1^2&\rho\sigma_1\sigma_2\\\rho\sigma_1\sigma_2&\sigma_2^2 \end{pmatrix}$. 
Then 
$$EX^2Y^2=\frac{\partial^4}{\partial^2t_1\,\partial^2t_2}\,M(t_1,t_2)\Big|_{t-1=t_2=0}=
4 \mu _1 \mu _2 \rho  \sigma _1 \sigma _2+\left(\mu _1^2+\sigma _1^2\right) \left(\mu
   _2^2+\sigma _2^2\right)+2 \rho ^2 \sigma _1^2 \sigma _2^2,
$$ 
whence 
$$Cov(X^2,Y^2)=4 \mu _1 \mu _2 \rho  \sigma _1 \sigma _2+2 \rho ^2 \sigma _1^2 \sigma _2^2. 
$$
Special cases of this formula are 
$$Var(X^2)=Cov(X^2,X^2)=4 \mu _1^2 \sigma _1^2+2 \sigma _1^4 
$$
and 
$$Var(Y^2)=Cov(Y^2,Y^2)=4 \mu _2^2 \sigma _2^2+2 \sigma _2^4.  
$$
So, the correlation between $X^2$ and $Y^2$ is 
$$\frac{Cov(X^2,Y^2)}{\sqrt{Var(X^2)Var(Y^2)}}
=\frac{4 \mu _1 \mu _2 \rho  \sigma _1 \sigma _2+2 \rho ^2 \sigma _1^2 \sigma
   _2^2}{\sqrt{\left(4 \mu _1^2 \sigma _1^2+2 \sigma _1^4\right) \left(4 \mu _2^2 \sigma
   _2^2+2 \sigma _2^4\right)}}
=\frac{2 \mu _1 \mu _2 \rho  \sigma _1 \sigma _2+ \rho ^2 \sigma _1^2 \sigma
   _2^2}{\sqrt{\left(2 \mu _1^2 \sigma _1^2+ \sigma _1^4\right) \left(2 \mu _2^2 \sigma
   _2^2+ \sigma _2^4\right)}}.
$$
In particular, the correlation between $X^2$ and $Y^2$ is $\rho^2$ when $\mu_1=\mu_2=0$.