The solid angle can be expressed in terms of complete elliptic integrals of first and third kind.
I believe this is first derived by F. Paxton in 1959.
THE REVIEW OF SCIENTIFIC INSTRUMENTS. VOLUME 30, NUMBER 4 APRIL, 1959.
Solid Angle Calculation for a Circular Disk. F. PAXTON
A copy of the article can be found
here.
Below is my attempt to derive the same result in an alternate manner.
The integral is indeed trickier than I thought.
To avoid confusion with the usage of $x,y,z$ as coordinates, let us change the problem and assume the circular disk $C$ we are dealing with lies in the plane $z = b$, centered at $(a,0,b)$ with radius $s$. We will assume $a, b > 0$ and $a \ne s$.
In terms of ordinary spherical polar coordinates $(r,\theta,\phi)$, points on the plane $z = b$ can be parametrized as:
$$(x,y,z) = (b\tan\theta \cos\phi,\,b\tan\theta\sin\phi,\, b)$$
Let $\rho = \sqrt{x^2+y^2} = b\tan\theta$, the "element" for integrating the solid angle is given by:
$$d\Omega
= \sin\theta d\theta \wedge d\phi
=\frac{\tan\theta d\tan\theta}{\sqrt{1 + \tan\theta^2}^3} \wedge d\phi
= \frac{b \rho d\rho \wedge d\phi}{\sqrt{\rho^2+b^2}^3}
= \frac{b\;dx \wedge dy}{\sqrt{\rho^2+b^2}^3}
$$
Introduce complex coorindates $\eta = x + iy$ and $\bar{\eta} = x - iy$, the solid angle
extended by $C$ can be rewritten as:
$$\begin{align}
\Omega_{C}
= & \int_{C} d\Omega
= \frac{b}{2i}\int_{C} \frac{d\bar{\eta} \wedge d\eta}{\sqrt{\eta\bar{\eta}+b^2}^3}
=\color{blue}{^{[1]}} \frac{b}{i} \int_C d\left(\frac{1}{\eta}\left( \frac{1}{b} - \frac{1}{\sqrt{\eta\bar{\eta}+b^2}}\right)\right) \wedge d\eta\\
= & \frac{b}{i} \int_{\partial C} \left(\frac{1}{b} - \frac{1}{\sqrt{\eta\bar{\eta} + b^2}}\right) \frac{d\eta}{\eta}
= 2\pi \delta_{C} + ib \int_{\partial C} \frac{d\eta}{\eta \sqrt{\eta\bar{\eta} + b^2}}
\end{align}$$
where $\delta_{C} = 1 \text{ or } 0$ depends on whether $s > a$ or $< a$. In other words, whether $\partial C$ contains $0$ or not.
On $\partial C$, we can parametrize $\eta$ as $s e^{it} + a$ and $\bar{\eta}$ as $s e^{-it} + a$. Substitute this into above expression, we obtain:
$$\begin{align}
&\Omega_C - 2\pi \delta_{C}\\
= &ib \int_{-\pi}^{\pi} \frac{ dt }{\sqrt{ s^2 + a^2 + b^2 + 2sa\cos t}}\frac{ is e^{it}}{s e^{it} + a }\\
= & -b
\int_{0}^{\pi} \frac{ dt }{\sqrt{ s^2 + a^2 + b^2 + 2sa\cos t}}
\left( \frac{s e^{it}}{s e^{it} + a } + \frac{s e^{-it}}{s e^{-it} + a } \right)\\
= & -b
\int_{0}^{\pi} \frac{ dt }{\sqrt{ s^2 + a^2 + b^2 + 2sa\cos t}}
\frac{2s(s + a\cos t)}{s^2 + a^2 + 2sa\cos t}\\
= & -b
\int_{0}^{\pi} \frac{ dt }{\sqrt{ s^2 + a^2 + b^2 + 2sa\cos t}}
\left( 1 + \frac{s^2 - a^2}{s^2 + a^2 + 2sa\cos t} \right)\\
= & -2b
\int_{0}^{\pi/2} \frac{ dt }{\sqrt{ s^2 + a^2 + b^2 + 2sa\cos(2t)}}
\left( 1 + \frac{s^2 - a^2}{s^2 + a^2 + 2sa\cos(2t)} \right)\\
= & -2b
\int_{0}^{\pi/2} \frac{ dt }{\sqrt{ (s+a)^2 + b^2 - 4sa\sin^2(t)}}
\left( 1 + \frac{s^2 - a^2}{(s+a)^2 - 4sa\sin^2(t)} \right)
\end{align}$$
The last integral can be expressed in terms of the
complete elliptic integral
of first and third kind:
$$
K(k) = \int_{0}^{\frac{\pi}{2}} \frac{dt}{\sqrt{1-k^2\sin^{2}t}}
\quad\text{ and }\quad
\Pi(n,k) = \int_{0}^{\frac{\pi}{2}} \frac{dt}{(1-n\sin^{2}t)\sqrt{1-k^2\sin^{2}t}}
$$
The final results are:
$$\Omega_{C} = 2\pi\delta_C -\frac{2b}{\sqrt{(s+a)^2+b^2}}
\left(\;K(k) + \left(\frac{s-a}{s+a}\right) \Pi(n,k)\;\right)
$$
where $\displaystyle \;\;k = \sqrt{\frac{4sa}{(s+a)^2+b^2}}
\quad\text{ and }\quad
n = \frac{4sa}{(s+a)^2}
$.
Notes
$\color{blue}{[1]}$ When $C$ contains $(0,0,b)$, we need the $\frac{1}{b}$ term. It regularize the 1-form at $\eta = 0$ and make Stroke's theorem continue to work.
Both of your methods are correct, and the flux through the sphere being $0$ is actually what we would expect, as your field is purely rotational and therefore the field vectors all point along the surface of the sphere. We can see this by observing:
$$\nabla \cdot \vec{F}=\frac{\partial(-y)}{\partial x}+\frac{\partial x}{\partial y}=0$$
And:
$$\nabla \times \vec{F}=\begin{vmatrix}\boldsymbol{\hat{\imath}} & \boldsymbol{\hat{\jmath}} & \boldsymbol{\hat{k}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ -y & x & 0\end{vmatrix}=2\boldsymbol{\hat{k}}$$
And so we can see that the field behaves purely rotationally, if we look at the vector plot of your vector field, this becomes more clear:
Best Answer
PRIMER:
Note that in regions that are absent of charge, $\nabla \cdot E=0$. Let $V$ be a region, bounded by the surface $S$, that contains a point charge. Furthermore, let $V'$ be a region bounded by the surface $S'$, contained in $V$.
If the point charge is in $V'$, then it is outside of the region $V\setminus V'$ and from Gauss's law and the Divergence Theorem we have
$$\begin{align} \int_{V\setminus V'}0\,dV&=\int_{V\setminus V'}\nabla\cdot\vec E\,dV\\\\ &=\oint_{S} \vec E\cdot \hat n\,dS-\oint_{S'} \vec E\cdot \hat n\,dS \end{align}$$
from which we conclude
$$\oint_{S} \vec E\cdot \hat n\,dS=\oint_{S'} \vec E\cdot \hat n\,dS$$
That is, the value of the surface integral is independent of the surface over which the integration is taken.
If the point charge is in $V\setminus V'$, then $\nabla \cdot E=0$ in $V'$ and $\oint_{S'}\vec E\cdot \hat n\,dS=0$.
We now will show explicitly that if $S'$ is any sphere that contains a point charge, that $\oint_{S'}\vec E\cdot \hat n\,dS=q/\varepsilon_0$ (We could have even chosen a sphere for which the point charge is at the origin.). To that end we proceed.
The electric field $\vec E(\vec r')$ due to a point charge at $\vec r$ at points on the surface of a sphere, centered at $0$ and with radius $a$ is
$$\vec E(\vec r')=\frac{q(a\hat r'-\vec r)}{4\pi\varepsilon_0|\vec r-a\hat r'|^3}$$
So, we wish to compute the integral
$$\oint_{|\vec r'|=a} \vec E(\vec r')\cdot \hat n'\,dS'=q\int_0^{2\pi}\int_0^\pi \left(\frac{a\hat r'-\vec r}{4\pi \epsilon_0|\vec r-a\hat r'|^3}\right)\cdot\hat r'\,a^2\,\sin(\theta')\,d\theta'\,d\phi'$$
Due to the spherical symmetry, without loss of generality, we align the polar axis so that the point charge lies at $\vec r=r\hat z$. Then,
$$\begin{align} \oint_{|\vec r'|=a} \vec E(\vec r')\cdot \hat n'\,dS'&=\frac{qa^2}{2\varepsilon_0}\int_0^\pi \frac{a-r\cos(\theta')}{(r^2+a^2-2ar\cos(\theta'))^{3/2}}\sin(\theta')\,d\theta'\\\\ &=\frac{qa^2}{2\varepsilon_0}\int_{-1}^1 \frac{a-rx}{(r^2+a^2-2arx)^{3/2}}\,dx\tag1 \end{align}$$
The integral on the right-hand side of $(1)$ is
$$\begin{align}\int_{-1}^1 \frac{a-rx}{(r^2+a^2-2arx)^{3/2}}\,dx&=-\left.\left(\frac{r-ax}{a^2(r^2+a^2-2arx)^{1/2}}\right)\right|_{-1}^1\\\\ &=\frac{r+a}{a^2(r+a)}-\frac{r-a}{a^2|r-a|}\\\\ &=\begin{cases} 0&, r>a\\\\ \frac 2{a^2}&,r<a\tag2 \end{cases} \end{align}$$
Using $(2)$ in $(1)$ we find the coveted result
$$\oint_{|\vec r'|=a} \vec E(\vec r')\cdot \hat n'\,dS'=\begin{cases}\frac{q}{\varepsilon_0}&,r<a\\\\ 0&,r>a\end{cases}$$
as was to be shown!