I have a circle in the $yz$-plane centered at $(0,2,0)$ with radius $1$.
The surface $\Sigma$ is obtained by rotating the circle around the $z$-axis.
I want calculate the flux of the vector field $F(x,y,z)=(\sqrt{y^2+z^2},\arctan(x+z),\cfrac{(z+1)^2}{\sqrt{x^2+y^2}})$
We can parametrize the circle with $\gamma(t)=(0,2+\sin{t},\cos{t})$ with $t \in [0,2\pi]$ so the surface $\Sigma=(-\sin{v}(2+\sin{u}), \cos{v}((2+\sin{u}), \cos{u})$ with $u,v \in [0,2\pi]$.
Can I use Stoke's theorem to solve this problem?
No you cannot use the Stoke's theorem for this problem. You can however use the divergence theorem for this type of problem though. It states that:
$$flux = \iint_\Sigma F\cdot dS = \iiint \nabla \cdot F dV$$
which essentially translates to the flux of the field $F$ through the closed surface is equal to the divergence of the field $F$ through the entire volume enclosed by that surface.
Note that the RHS is easier to calculate than the LHS because
$$\begin{aligned} \cfrac{\partial F_1}{\partial x} &= 0\\ \cfrac{\partial F_2}{\partial y} &= 0\\ \cfrac{\partial F_3}{\partial z} &= \cfrac{2(z+1)}{\sqrt{x^2+y^2}} \end{aligned}$$
Now, we have to parameterize the volume of the torus (not the area as you have already done). This is easily done by further noting that parameterizing the area of circle means that we have a radius $a \in [0,1]$ as opposed to $a=1$ on the boundary.
Thus, $\gamma(u,a) = (0, 2+a\sin{u},a\cos{u})$ with $u \in [0,2\pi), a \in [0,1]$. And the volume $V(u,a,v)=(-\sin{v}(2+a \sin{u}), \cos{v}(2+a \sin{u}), a \cos{u})$ with $u \in [0,2\pi), v \in [0,2\pi), a \in [0,1]$.
EDIT: Before we revisit the integral, we first need to compute the Jacobian for this transformation from $(x,y,z)$ coordinates to $(u,v,a)$ coordinates.
$$\begin{aligned} Jacobian&=\begin{vmatrix}\cfrac{\partial x}{\partial u} & \cfrac{\partial x}{\partial v} & \cfrac{\partial x}{\partial a}\\ \cfrac{\partial y}{\partial u} & \cfrac{\partial y}{\partial v} & \cfrac{\partial y}{\partial a}\\ \cfrac{\partial z}{\partial u} & \cfrac{\partial z}{\partial v} & \cfrac{\partial z}{\partial a} \end{vmatrix}\\ &= \begin{vmatrix}-a\cos{u}\sin{v} & -\cos{v}(2+a \sin{u}) & -\sin{u}\sin{v}\\ a \cos{u} \cos{v} & -\sin{v}(2+a \sin{u}) & \sin{u}\cos{v}\\ -a\sin{u} & 0 & \cos{u} \end{vmatrix}\\ &= a\cos^2{u}\sin^2{v}(2+a\sin{u}) + a\sin^2{u} \cos^2{v}(2+a\sin{u})+a \sin^2{u}\sin^2{v}(2+a\sin{u})+a\cos^2{u}\cos^2{v}(2+a\sin{u})\\ &= a\sin^2{v}(2+a\sin{u})+a\cos^2{v}(2+a\sin{u})\\ &= a(2+a\sin{u}) \end{aligned}$$
Then we can rewrite our volume differential using this Jacobian as $dxdydz=a(2+a\sin{u})dudadv$.
$$\begin{aligned} flux&=\iiint \cfrac{2(z+1)}{\sqrt{x^2+y^2}} dV\\ &=\int_0^{2\pi}\int_0^1 \int_0^{2\pi} \cfrac{2a \cos{u} + 2}{a\sin{u}+2}(a(2+a\sin{u}))du da dv \\ &=\int_0^{2\pi}\int_0^1 \int_0^{2\pi} (2a^2 \cos{u} + 2a)du da dv \\ &=\int_0^{2\pi}\int_0^1 4 \pi a da dv\\ &=\int_0^{2\pi} 2 \pi dv\\ &=4\pi^2 \end{aligned}$$
Note how much easier the RHS of the divergence theorem is compared to the monstrosity that is the LHS in this case.