Double Integral $\int\limits_0^a\int\limits_0^a\frac{dx\,dy}{(x^2+y^2+a^2)^\frac32}$

Both the function that you are integrating as the region over which you are integrating it get unchanged if you exchange $x$ with $y$. Therefore, your integral is equal to$$2\int_0^a\int_0^x\frac1{(a^2+x^2+y^2)^{3/2}}\,\mathrm dy\,\mathrm dx.$$You can compute this integral using polar coordinates: $\theta$ can take values in $\left[0,\frac\pi4\right]$ and, for each $\theta$, $r$ can take values in $\left[0,\frac a{\cos\theta}\right]$. And\begin{align}\int_0^{\pi/4}\int_0^{a/\cos(\theta)}\frac r{(a^2+r^2)^{3/2}}\,\mathrm dr\,\mathrm d\theta&=\int_0^{\pi/4}\frac{1-\frac1{\sqrt{\sec ^2(\theta)+1}}}a\,\mathrm d\theta\\&=\frac1a\left(\frac\pi4-\int_0^{\pi/4}\frac1{\sqrt{\sec ^2(\theta)+1}}\,\mathrm d\theta\right)\\&=\frac\pi{12a}.\end{align}Note that the final equality is equivalent to$$\int_0^{\pi/4}\frac1{\sqrt{\sec^2(\theta)+1}}\,\mathrm d\theta=\frac\pi6.$$This can be justified as follows: you do $\theta=\arccos\left(\sqrt x\right)$ and $\mathrm d\theta=-\frac1{2\sqrt{x-x^2}}\,\mathrm dx$. Doing this, you will get\begin{align}\int_0^{\pi/4}\frac1{\sqrt{\sec^2(\theta)+1}}\,\mathrm d\theta&=\int_1^{1/2}-\frac1{2\sqrt{1-x^2}}\\&=\frac12\int_{1/2}^1\frac1{\sqrt{1-x^2}}\,\mathrm dx\\&=\frac12\left(\arcsin\left(1\right)-\arcsin\left(\frac12\right)\right)\\&=\frac12\left(\frac\pi2-\frac\pi6\right)\\&=\frac\pi6.\end{align}

Let $y=\sqrt{a^2+x^2} \tan \theta$, $dy=\sqrt{a^2+x^2}~ \mbox{sec}^2 \theta~ d\theta$ Then $$\int_{0}^{a} \int_{0}^{a} \frac{dx ~ dy} {(a^2+x^2+y^2)^{3/2}} =\int_{0}^{a}\frac{dx}{a^2+x^2}\int_{0}^{\tan^{-1}(a/\sqrt{a^2+x^2})} \cos \theta ~d\theta$$ $$ = a \int_{0}^{a} \frac{dx}{(x^2+a^2)\sqrt{(2a^2+x^2)}}.$$ $$ \Rightarrow I=\int_{0}^{a}\frac{1}{2i} \frac{dx}{\sqrt{2a^2+x^2}} \left [ \frac{1}{x-ia}-\frac{1}{x+ia}\right]=\Im \left( \int_{0}^{a}\frac{dx}{(x-ia)\sqrt{2a^2+x^2}}\right).$$ Use $(x-ia)=1/t$, then $$I=-\Im \left (\int_{i/a}^{(1+i)/(2a)} \frac{dt}{\sqrt{a^2t^2+2iat+1}}\right) =-\Im \int_{2i}^{(1+3i)/2} \left(\frac{dv}{a\sqrt{v^2+(\sqrt{2})^2}}\right).$$

Letting $v=\sqrt{2} \tan \phi$, we find that $$I=-\Im \left( a^{-1} \ln \left[ \frac{1+3i+\sqrt{6i}}{2i \sqrt{2}(1+\sqrt{2})}\right]\right)= \frac{\pi}{6a}.$$

The natural course of action whenever you see $x^2 + y^2$ in a multiple integration problem is to convert to polar coordinates.

Set $x = r\cos\theta$, $y = r\sin\theta$, $0 \leq \theta <2\pi$. Then $dx~dy = r~dr~d\theta$, and the integrand becomes $$ \frac{r}{(r^2 + a^2)^{3/2}}~dr~d\theta. $$ This integral doesn't depend on $\theta$, so in fact when you integrate in $r$ it is a one-variable integral, and from single-variable calculus we know that the substitution $u = r^2 + a^2$ will do the trick. The only thing left to do is to find the limits of integration.

The region $0 \leq x \leq a$, $0 \leq y \leq a$ is a square of side length $a$ in the first quadrant of the plane. This can be parametrized by $0 \leq \theta \leq \frac{\pi}{4}$, $0 \leq r \leq \sqrt{a^2 + (a\sin\theta)^2} = a\sqrt{1 + \sin^2 \theta}$. (Draw a picture! This is the meat of the problem.) Therefore the integral becomes $$ \int_{\theta = 0}^{\pi/4}\int_{r=0}^{a\sqrt{1+\sin^2\theta}} \frac{r}{(r^2 + a^2)^{3/2}}~dr~d\theta. $$

Can you take it from here?