Does $\iint_D \frac{x^2}{x^2+y^2} dx dy $ converge on $D= \left\{ (x, y) : x^2+y^2\leq ax \right\} $ ? If yes, what value does it converge to? Powers equal to $2$ entice the polar substitution.
$D= \left\{ (x, y) : x^2+y^2\leq ax \right\}$, so $0 \leq r \leq a \cos \phi.$ For the domain to make sense, we need either $\phi \in [0, \frac{\pi}{2}] \cup [\frac{3 \pi}{2}, 2\pi)$ and $a \geq 0$, or $\phi \in [\frac{\pi}{2}, \frac{3 \pi}{2}]$ and $a \leq 0$.
Everything below was achieved with help from the comments.
$$\iint_D \frac{x^2}{x^2+y^2} dx dy = \iint_{D'} r \cos^2 \phi \ d\phi dr.$$
Case $1$:
Since $(0, 0) \in D $ (equivalently, $r=0$ for any $\phi$) makes the integrand indefinite, 
$$
\lim_{\epsilon \rightarrow 0}
\int_{0}^{\frac{\pi}{2}} \cos^2 \phi \ d \phi 
\int_{\epsilon}^{a \cos \phi} r \ dr 
+ 
\lim_{\epsilon \rightarrow 0}
\int_{\frac{3 \pi}{2}}^{2\pi} \cos^2 \phi \ d \phi 
\int_{\epsilon}^{a \cos \phi} r \ dr
$$
$$
2\lim_{\epsilon \rightarrow 0}
\int_{0}^{\frac{\pi}{2}} \cos^2 \phi \ d \phi 
\int_{\epsilon}^{a \cos \phi} r \ dr 
$$
$$
2\lim_{\epsilon \rightarrow 0}
\int_{0}^{\frac{\pi}{2}} \left( \frac{a^2 \cos^2 \phi}{2}-\frac{\epsilon^2}{2} \right) \cos^2 \phi \ d \phi 
$$
$$
\lim_{\epsilon \rightarrow 0}
\int_{0}^{\frac{\pi}{2}} \left( {a^2 \cos^2 \phi}-{\epsilon^2} \right) \cos^2 \phi \ d \phi 
$$
Put $0$ instead of $\epsilon$, because it won't let me integrate:
$$
a^2
\int_{0}^{\frac{\pi}{2}} { \cos^4 \phi} \ d \phi = \frac{a^2}{2} Β \left(\frac{1}{2}, \frac{5}{2} \right) = \frac{a^2}{2} \cdot \frac{\Gamma(\frac{1}{2}) \Gamma(\frac{5}{2})}{\Gamma(3)} = \frac{a^2}{2} \cdot \frac{\sqrt{\pi} \cdot 0.75 \sqrt{\pi}}{2} = \frac{3 a^2 \pi}{16}
$$
For case $2$, the result should be same.
Was this correct?
Formula $Β(m, n)=2\int_0^{\pi/2}\sin^{2m-1} \theta \cos^{2n-1} \theta \ d \theta$ was used, but I'm not sure I can extract it myself, and would be grateful for reference. 
The answer would be, that for fixed $a$, it's convergent. (I am not sure about the definition of convergence here, and can't look it up, since the classes were missed, and google is blocked). 

For Oliver Jones:
let $r_i=a\frac{i}{i+1}$, $\phi_i=\frac{\delta}{i}.$
$$\iint_{D_\delta} r \cos^2 \phi \ d\phi dr= \lim_{i \rightarrow \infty} \sum_i \left(a\frac{i}{i+1}\right) \left(\cos^2 \frac{\delta}{i} \right) 
\left( 
a\frac{i+1}{i+2} - a\frac{i}{i+2}
\right) 
\left( 
\frac{\delta}{i+1}-\frac{\delta}{i}
\right) =
\lim_{i \rightarrow \infty} \sum_i \left(a\frac{i}{i+1}\right) \left(1 \right) 
\left( 
\frac{a}{(i+1)(i+2)}
\right) 
\left( 
\frac{-\delta}{i(i+1)}
\right)
=
 \sum_{i=0}^{\infty} \left(a^2\frac{1}{i+1}\right) \left(1 \right) 
\left( 
\frac{1}{(i+1)(i+2)}
\right) 
\left( 
\frac{-\delta}{(i+1)}
\right).
$$
It converges because of the $4$th power in the denominator.
Or maybe like this:
$$
\lim_{\epsilon \rightarrow 0} \int_{\epsilon}^a x \ dx \int_{-\sqrt{x(a-x)}}^{\sqrt{x(a-x)}} \frac{1}{1+ (y/x)^2} d (y/x)$$
$$
\lim_{\epsilon \rightarrow 0} \int_{\epsilon}^a 2 x \arctan \sqrt{x(a-x)} \ dx  $$
How to proceed here, I am not sure.

$$
\lim_{\delta \rightarrow 0}\iint_{\delta^2 \leq r^2+ar \cos \theta + a^2/4 \leq a^2/4} r \cos ^2 \theta \ dr \ d\theta
$$
Let's choose $\theta \in [-\pi/2, \pi/2], a\geq 0$. 
Then 
$
\delta^2 \leq r^2 - ar + a^2/4 \leq  r^2+ar \cos \theta + a^2/4 \leq r^2 + a^2/4 \leq a^2/4.
$
$
0 \leq r^2 - ar + a^2/4 -\delta^2 \Longrightarrow r\in (0, \frac{a-2\delta}{2}).
$
$$
\lim_{\delta \to 0} \int_{-\pi/2}^{\pi/2} \cos ^2 \theta \ d \theta \int_0^{(a-2\delta)/2} r dr = \frac{a^2}{8} В(1/2, 3/2) = \frac{a^2 \pi}{16}
$$

I get the same result as your first one. But I don't see the need for special functions, or introducing the $\epsilon$ you use.
Consider the region $x^2+y^2\leq ax$. This is the same as $$\left(x-\frac{a}{2}\right)^2+y^2\leq\frac{a^2}{4}$$
So the region is a disc centered at $\left(\frac{a}{2},0\right)$ with its boundary touching the origin. I'll assume $a>0$. A disc like this is described in polar coordinates by $r\leq a\cos(\theta)$ for $\theta\in[-\pi/2,\pi/2]$.
So the integral is 
$$
\begin{align}
&\int_{\theta=-\pi/2}^{\theta=\pi/2}\int_{r=0}^{r=a\cos(\theta)}r\cos^2(\theta)\,dr\,d\theta\\
&=\int_{\theta=-\pi/2}^{\theta=\pi/2}\cos^2(\theta)\int_{r=0}^{r=a\cos(\theta)}r\,dr\,d\theta\\
&=\int_{\theta=-\pi/2}^{\theta=\pi/2}\cos^2(\theta)\frac{1}{2}a^2\cos^2(\theta)\,d\theta\\
&=\frac{a^2}{2}\int_{\theta=-\pi/2}^{\theta=\pi/2}\cos^4(\theta)\,d\theta\\
&=\frac{a^2}{2}\int_{\theta=-\pi/2}^{\theta=\pi/2}\left(\frac{1+\cos(2\theta)}{2}\right)^2\,d\theta\\
&=\frac{a^2}{4}\int_{\varphi=-\pi}^{\varphi=\pi}\left(\frac{1+\cos(\varphi)}{2}\right)^2\,d\varphi\\
&=\frac{a^2}{4}\int_{\varphi=-\pi}^{\varphi=\pi}\left(\frac14+\cos(\varphi)+\frac{1}{4}\cos^2(\varphi)\right)\,d\varphi\\
&=\frac{a^2}{4}\int_{\varphi=-\pi}^{\varphi=\pi}\left(\frac14+\frac{1}{4}\cos^2(\varphi)\right)\,d\varphi\\
&=\frac{a^2}{4}\int_{\varphi=-\pi}^{\varphi=\pi}\left(\frac14+\frac{1}{8}\left(1+\cos(2\varphi)\right)\right)\,d\varphi\\
&=\frac{a^2}{4}\int_{\varphi=-\pi}^{\varphi=\pi}\left(\frac38+\frac{1}{8}\cos(2\varphi)\right)\,d\varphi\\
&=\frac{a^2}{4}\int_{\varphi=-\pi}^{\varphi=\pi}\frac38\,d\varphi\\
&=\frac{3a^2}{32}(2\pi)=\frac{3a^2\pi}{16}
\end{align}$$
