Calculating a flux integral 
Let $$F=(xe^{xy}-2xz+2xy\cos^2 z, y^2\sin^2 z-y e^{xy}+y, x^2+y^2+z^2)$$ and $V$ be the solid in space bounded by $z=9-x^2-y^2$ and $z=0$. I am trying to compute the flux integral $\iint_{\partial V}F\cdot n \ dS$, $n$ being the outward unit normal.

Setting $r(x,y)=(x,y,9-x^2-y^2)$, I found that $r_x\times r_v=(2x,2y,1)$ and $$\iint_{\partial V}F\cdot n \ dS=\iint_D F\cdot r_x\times r_y \ dA$$ where $$F\cdot r_x\times r_y=2e^{xy}(x^2-y^2)-36x^2+4x^4+4x^2y^2+4x^2y\cos^2(9-x^2-y^2)+2y^3\sin^2(9-x^2-y^2)+2y^2;$$ after the substitution $x=r\cos t, \ y=r\sin t$ I get $$f(r,t)=r(F\cdot r_x\times r)=2r^3e^{r^2\sin t \cos t}(\cos^2 t-\sin^2 t)-36r^3\cos^2 4r^5\cos^2 t + \\4r^4\cos^2 t\sin t \cos^2{(9-r^2)}+2r^4\sin^3 t \sin^2(9-r^2)+2r^3\sin^2 t $$ and I need to compute $$\int_0^{2\pi}\int_0^3 f(r,t)drdt$$
I wonder whether I can further simplify $f(r,t)$? The current expression looks to cumbersome and it seems like a hassle to compute the integral if no simplifications can be made.

Why aren't you using the divergence theorem? Compute $$\begin{align} {\rm div}\,F(x,y,z) &= \partial_x(xe^{xy}-2xz+2xy\cos^2 z)+\partial_y(y^2\sin^2 z-y e^{xy}+y)+\partial_z(x^2+y^2+z^2) \\ &= \require{cancel} \cancel{e^{xy}} + \cancel{xye^{xy}} -
 \cancel{2z} + 2y\cos^2z + 2y\sin^2 z - \cancel{e^{xy}} - \cancel{xye^{xy}}+1+\cancel{2z} \\ &= 2y + 1.\end{align}$$Then $$\begin{align} \iint_{\partial V} F\cdot n\,{\rm d}S &= \iiint_V {\rm div}\,F(x,y,z)\,{\rm d}x\,{\rm d}y\,{\rm d}z \\ &= \iiint_V 1+2y \,{\rm d}x\,{\rm d}y\,{\rm d}z \\ &= \int_0^{2\pi}\int_0^3\int_0^{9 - r^2} (1+2r\sin \theta)r\,{\rm d}z\,{\rm d}r\,{\rm d}\theta, \end{align}$$which is easy. Can you finish?
