Exact expression of a trigonometric integral

With CAS help:

$$\int _0^{\pi }\int _0^{\pi }\frac{\sin ^2(x) \sin ^2(y)}{a+\cos (x)+\cos (y)}dydx=\\\mathcal{L}_q\left[\int _0^{\pi }\int _0^{\pi }\mathcal{L}_a^{-1}\left[\frac{\sin ^2(x) \sin ^2(y)}{a+\cos (x)+\cos (y)}\right](q)dydx\right](a)=\\\mathcal{L}_q\left[\int_0^{\pi } \left(\int_0^{\pi } e^{-q (\cos (x)+\cos (y))} \sin ^2(x) \sin ^2(y) \, dx\right) \, dy\right](a)=\\\mathcal{L}_q\left[\int_0^{\pi } \frac{e^{-q \cos (y)} \pi I_1(q) \sin ^2(y)}{q} \, dy\right](a)=\\\mathcal{L}_q\left[\frac{\pi ^2 I_1(q){}^2}{q^2}\right](a)=\\\frac{a \pi ^2}{2}-\frac{2}{3} a \pi E\left(\frac{4}{a^2}\right)-\frac{1}{6} a^3 \pi E\left(\frac{4}{a^2}\right)-\frac{2}{3} a \pi K\left(\frac{4}{a^2}\right)+\frac{1}{6} a^3 \pi K\left(\frac{4}{a^2}\right)=\\\frac{\pi ^2 \, _3F_2\left(\frac{1}{2},1,\frac{3}{2};2,3;\frac{4}{a^2}\right)}{4 a}$$

for: $a>2$

$$\frac{1}{6} \pi \left(-2 \left(a^2-4\right) K\left(\frac{a^2}{4}\right)-2 \left(a^2+4\right) E\left(\frac{a^2}{4}\right)+3 \pi a\right)$$

for: $a<2$

where: $K$,$E$ gives the elliptic integral of the first kind and second kind.

Mathematica code:

HoldForm[Integrate[(Sin[x]^2*Sin[y]^2)/(a + Cos[x] + Cos[y]), {x, 0, Pi}, {y, 0, Pi}] == (a \[Pi]^2)/2 - 2/3 a \[Pi] EllipticE[4/a^2] - 1/6 a^3 \[Pi] EllipticE[4/a^2] - 2/3 a \[Pi] EllipticK[4/a^2] + 1/6 a^3 \[Pi] EllipticK[4/a^2] == Pi^2/(4 a)*HypergeometricPFQ[{1/2, 1, 3/2}, {2, 3}, 4/a^2]] // TeXForm

Plot a solution:

f[a_?NumericQ] := NIntegrate[(Sin[x]^2*Sin[y]^2)/(a + Cos[x] + Cos[y]), {x, 0, Pi}, {y,0, Pi}]; 
g[a_] := (a \[Pi]^2)/2 - 2/3 a \[Pi] EllipticE[4/a^2] - 
1/6 a^3 \[Pi] EllipticE[4/a^2] - 2/3 a \[Pi] EllipticK[4/a^2] + 
1/6 a^3 \[Pi] EllipticK[4/a^2]; Plot[{f[a], g[a]}, {a, 2, 20}, 
PlotStyle -> {Red, {Dashed, Black}}, PlotLabels -> {"integral", "Analytic solution"}]

Partial answer:

As the integrand is an even function, write $$I(a)=\frac14\int_{-\pi}^\pi\int_{-\pi}^\pi\frac{\sin^2x\sin^2y}{a+\cos x+\cos y}\,dx\,dy$$ and let $w=e^{ix}$. Then for $C_1=\{w:|w|=1\}$, we have \begin{align}\int_{-\pi}^\pi\frac{\sin^2x}{a+\cos x+\cos y}\,dx&=\oint_{C_1}\frac{\frac{(w-1/w)^2}{-4}}{b+\frac{w+1/w}2}\frac{dw}{iw}=\frac i2\oint_{C_1}\frac{(w^2-1)^2}{w^2(w^2+2bw+1)}\,dw\end{align} where $b:=a+\cos y$. Inside $C_1$, there is a double pole at $0$, so letting $f(w)$ be the integrand, \begin{align}\operatorname{Res}(f,0)=\lim_{w\to0}\frac d{dw}(w^2f(w))=-2b.\end{align} Solving the quadratic in the denominator of $f$ gives $w=-b\pm\sqrt{b^2-1}$ but only the positive root lies in $C_1$. As this is a simple pole, \begin{align}\operatorname{Res}(f,-b+\sqrt{b^2-1})=\lim_{w\to-b+\sqrt{b^2-1}}(w+b-\sqrt{b^2-1})f(w)=2\sqrt{b^2-1}.\end{align} Hence \begin{align}\int_{-\pi}^\pi\frac{\sin^2x}{a+\cos x+\cos y}\,dx=\frac i2\cdot2\pi i(2\sqrt{b^2-1}-2b)=-2\pi(\sqrt{b^2-1}-b)\end{align} so that \begin{align}I(a)&=-\frac\pi2\int_{-\pi}^\pi\sin^2y\left(\sqrt{(a+\cos y)^2-1}-a-\cos y\right)\,dy\\&=-\frac\pi2\int_{-\pi}^\pi\sin^2y\sqrt{(a+\cos y)^2-1}\,dy+\frac{a\pi^2}2.\end{align} Now let $z=e^{iy}$. Then for $C_2=\{z:|z|=1\}$, we have \begin{align}\int_{-\pi}^\pi\sin^2y\sqrt{(a+\cos y)^2-1}\,dy&=\oint_{C_2}\left(\frac{z-1/z}{2i}\right)^2\sqrt{\left(a+\frac{z+1/z}2\right)^2-1}\cdot\frac{dz}{iz}\\&=\frac i8\oint_{C_2}\frac{(z^2-1)^2}{z^4}\sqrt{(z^2+2az+1)^2-4z^2}\,dz.\end{align} Evidently, there is a pole of order $4$ at $0$ which lies in $C_2$. However, due to the presence of $\sqrt\cdot$, branch points also need to be considered, and solving the equation $z^2+2az+1=\pm2z$ reveals that only $z=\pm1-a+\sqrt{a^2\mp2a}$ are the two branch points in $C_2$.

We use the identities $$ 2\sin \phi \sin \psi=\cos(\phi-\psi)-\cos(\phi+\psi) $$ and $$ \cos \phi+\cos \psi=2\cos\left(\frac{\phi+\psi}{2}\right)\cos\left(\frac{\phi-\psi}{2}\right) $$ to get $$ I=\int^{\pi}_{0}\int^{\pi}_{0}\frac{\sin^2 x \sin^2 y}{a+\cos x+\cos y}dxdy= $$ $$ \frac{1}{4}\int^{\pi}_{0}\int^{\pi}_{0}\frac{(\cos(x-y)-\cos(x+y))^2}{a+2\cos\left(\frac{x+y}{2}\right)\cos\left(\frac{x-y}{2}\right)}dxdy $$ Now we use $\cos \phi=2\cos^2(\phi/2)-1$ and arrive to $$ I=\int^{\pi}_{0}\int^{\pi}_{0}\frac{(\cos^2\left(\frac{x-y}{2}\right)-\cos^2\left(\frac{x+y}{2}\right))^2}{a+2\cos\left(\frac{x+y}{2}\right)\cos\left(\frac{x-y}{2}\right)}dxdy= $$ $$ =a^{-1}\int^{\pi}_{0}\int^{\pi}_{0}\frac{(\cos^2 A-\cos^2 B)^2}{1+2/a \cos A\cos B}dAdB. $$ Then set $t=\cos A$, $s=\cos B$ and use $$ \frac{d\left(\cos^{(-1)}(x)\right)}{dx}=-\frac{1}{\sqrt{1-x^2}} $$ to get $$ I=1/a\int^{1}_{-1}\int^{1}_{-1}\frac{(t^2-s^2)^2}{1+(2/a)ts}\frac{1}{\sqrt{1-t^2}}\frac{1}{\sqrt{1-s^2}}dtds $$ If $|2/a|\leq 1$ and $|t|,|s|<1$, we expand $(1+(2/a)ts)^{-1}$ into $\sum^{\infty}_{k=0}(-1)^k(2/a)^kt^ks^k$ we can split the two integrals (that of $t$ and that of $s$). We set $$ I_k:=1/a\int^{1}_{-1}\int^{1}_{-1}(t^2-s^2)^2(-1)^k(2/a)^kt^ks^k\frac{1}{\sqrt{1-t^2}}\frac{1}{\sqrt{1-s^2}}dtds. $$ Then $I_{2k+1}=0$ and $$ I_{2k}=\frac{2^{2k}(1/a)^{2k+1}\pi\Gamma\left(\frac{2k+1}{2}\right)\Gamma\left(\frac{2k+3}{2}\right)}{\Gamma\left(k+2\right)\Gamma(k+3)}. $$ Hence summing $$ I=I(a)=\sum^{\infty}_{k=0}I_{2k}=\frac{a\pi^2}{2}-\frac{a\pi^2}{2}{}_2F_1\left(\frac{-1}{2},\frac{1}{2};2;\frac{4}{a^2}\right)= $$ $$ =\frac{a\pi^2}{2}-\frac{a(a^2+4)\pi}{6}E\left(4a^{-2}\right)+\frac{a(a^2-4)\pi}{6}K(4a^{-2}), $$ where $K(x)$,$E(x)$ are the complete elliptic integrals of the first and second kind. Examples of evaluations can be given if we consider the singular modulus $k_r$. For example with $a=2\sqrt{2}$, we get $k_1=1/2$ and we have $$ I\left(2\sqrt{2}\right)=\sqrt{2}\pi^2-2\sqrt{2\pi}\Gamma\left(\frac{3}{4}\right)^2-\frac{8}{3}\sqrt{2\pi}\Gamma\left(\frac{5}{4}\right)^2 $$ ...etc

see elliptic integrals singular valus