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.

If you are calculating the flux directly, then remember to normalize the normal to the surface. That is 
$$\mathbf{n}=\frac{\mathbf{r}_x {\times} \mathbf{r}_y}{\Vert {\mathbf{r}_x {\times} \mathbf{r}_y} \Vert}$$
HINT: A simpler approach is to use the Divergence theorem:
$$\iint_R{\mathbf{F}.\mathbf{n}}dS = \iiint_V{\nabla \cdot\mathbf{F}dV}$$
