I encountered this integral when trying to obtain a Fourier series for the function inside (in connection to this question).
Mathematica gives the following general solution (only valid for $|a|>1$), containing regularized generalized hypergeometric function:
$$A_n=\frac{2}{\pi} \int^{\pi}_0 \frac{\sin^2 (y)}{a+\cos(y)} \cos(ny) dy=$$
$$=\frac{4}{a-1} \left(\, _3\tilde{F}_2\left(1,\frac{3}{2},2;2-n, n+2;-\frac{2}{a-1}\right)-3 \, _3\tilde{F}_2\left(1,\frac{5}{2},3;3-n, n+3;-\frac{2}{a-1}\right)\right)$$
Some information about this particular function is here. I can't understand if the formulas on these page are general or not, and if I can apply them in my case.
However, I expect this expression to become simpler for the integer values of $n$. And indeed Mathematica give the explicit expressions up to $n=10$. The first few are:
$$A_0=2(a-\sqrt{a^2-1})$$
$$A_1=1-2a(a-\sqrt{a^2-1})$$
$$A_2=2 \sqrt{a^2-1}~~ (1-2a(a-\sqrt{a^2-1}))$$
$$A_3=2 \sqrt{a^2-1} ~~ ((4a^2-1)(a-\sqrt{a^2-1})-2a)$$
Technically, these $10$ expressions are enough for practical puproses, especially for the case $a>2$, because they become very small.
However, it would be great if there is some general simplification (so I can get rid of the special function).
Extend the integration domain to $[-\pi,\pi]$ using parity, then again using parity replace $\cos ny$ by $e^{in y}$, and finally make the change of variables $e^{iy}=z$. The result is an integral along the unit circle $|z|=1$ oriented counterclockwise: $$A_n=-\frac{1}{2\pi i}\oint_{|z|=1}\frac{\left(z-z^{-1}\right)^2 z^{n}dz}{\left(z+a+\sqrt{a^2-1}\right)\left(z+a-\sqrt{a^2-1}\right)}.$$ Now by residues we get $$A_{n\ge2}=-\operatorname{res}_{z=-a+\sqrt{a^2-1}}\frac{\left(z-z^{-1}\right)^2 z^{n}}{z+a+\sqrt{a^2-1}}=2(-1)^{n+1}\sqrt{a^2-1}\left(a-\sqrt{a^2-1}\right)^{n}.$$