(* We compute the flux of F through the paraboloid z=x^2+y^2 below the plane z=1 oriented down *)

F={x,y,z^4}

(* parametrization using cylindrical coordinates*)
q={r Cos[t],r Sin[t],r^2}
qr=D[q,r]
qt=D[q,t]

(* normal vector *)
n=Simplify[Cross[qt,qr]]

(* Double integral *)
Fcompq=F/.{x->q[[1]],y->q[[2]],z->q[[3]]}
g=Dot[Fcompq,n]
Integrate[Integrate[g,{r,0,1 }],{t,0,2Pi}]

(* alternate simpler parametrization, followed by change of variable *)
q={x,y ,x^2+y^2}
qx=D[q,x]
qy=D[q,y]

(* Normal vector *)
n=Simplify[Cross[qy,qx]]

Fcompq=F/.{x->q[[1]],y->q[[2]],z->q[[3]]}
g=Dot[Fcompq,n]

(* switch to polar coordinates to handle the region of integration (disk) *)
gg=Simplify[g/.{x ->r Cos[t],y->r Sin[t]}]
Integrate[Integrate[gg r,{r,0,1 }],{t,0,2 Pi}]