How to find outward-pointing normal vector for surface flux problems? Example problem included. Example Problem
Use a parameterisation to find the flux of $F = (3xy, 0, -z)$ outward (normal away from the z-axis) through the cone $z = 9\sqrt{x^2 + y^2}$, $0 \le z \le 9$.

My Work
We will use cylindrical coordinates to parameterise the cone.
The formula for the flux through a surface is $\iint_S F(S(\rho, \theta)) \cdot \overrightarrow{N} dt$ where $S(\rho, \theta)$ is the parameterised surface and $\overrightarrow{N} =  \dfrac{\partial{S}}{\partial{\rho}} \times \dfrac{\partial{S}}{\partial{\theta}}$.
The parameterised surface is $S(\rho, \theta) = (\rho \cos(\theta), \rho \sin(\theta), 9\rho) \forall \rho \in [0, 1], \theta \in [0, 2\pi]$.
$\dfrac{\partial{S}}{\partial{\rho}} = (\cos(\theta), \sin(\theta), 9)$
$\dfrac{\partial{S}}{\partial{\theta}} = (-\rho \sin(\theta), \rho \cos(\theta), 0)$
$\overrightarrow{N} = \dfrac{\partial{S}}{\partial{\rho}} \times \dfrac{\partial{S}}{\partial{\theta}} = (-9\rho \cos(\theta), -9\rho \sin(\theta), \rho)$
$F(S(\rho, \theta)) = (3\rho^2 \cos(\theta)\sin(\theta), 0, -9\rho)$
Flux $= \int^{2\pi}_{\theta = 0} \int^{1}_{0} (3\rho^2 \cos(\theta)\sin(\theta), 0, -9\rho) \cdot (-9\rho \cos(\theta), -9\rho \sin(\theta), \rho) d\rho d\theta$
$= -6\pi$

Question and Explanation 
I mostly understand flux problems, but one component that I have no idea how to do is to find the outward-pointing normal vector (as opposed to the inward-pointing normal vector). At the moment, I just take the cross product of the partial derivatives of the parameterised surface and hope I get the correct normal vector, which is obviously a serious deficiency in my knowledge.
How does one ensure that they get the outward-pointing normal vector, $\overrightarrow{N}$? I have no idea how to do this generally, across all (closed) surface flux problems. I would prefer an analytical method for finding this rather than graphing the surfaces every time, since I really don't think graphing everything is a robust/generalisable strategy -- especially in higher dimensions.
The above problem is an example of how I solve these problems. As you can see, I just calculated $\overrightarrow{N}$ with no way of knowing whether it was the outward-pointing normal vector. In this specific example, I have a feeling that it is actually the incorrect (inward-pointing) normal vector. This is an example of the deficiency in my skills.
I would greatly appreciate it if people could please take the time to explain a methodology for ensuring that I find the outward-pointing normal vector when solving these types of problems. I've been doing a lot of research on this, but I cannot find a clear explanation of a methodology. I would very much appreciate help from MSE members.
 A: It is a fact of life that $S^0$ consists of two points. This then implies that the normal line $n$ through a point $p$ of a hypersurface $S\subset{\mathbb R}^n$ has a natural origin, namely $p$, but two unit vectors spanning $n$, which differ by a factor $-1$. In a concrete geometrical or physical situation involving a two-dimensional surface $S\subset{\mathbb R}^3$ one has a particular "positive" normal ${\bf n}$ in mind, which one then describes as "outward", or "upward", etc., depending on the circumstances.
Assume that you are given a surface $S$ in geometrical or colloquial terms, say "a sphere with center $M$ and radius $r$, oriented outwards", and you then consult a catalogue of surface representations in order to obtain a parametric representation of $S$, then there is only a $50\%$ chance that the two parameters $u_1$, $u_2$ used generate via ${\bf f}_{.1}\times{\bf f}_{.2}$ the desired "orientation" of ${\bf n}$. For instance, in conection with spherical coordinates $(\phi,\theta)$ it plays a rôle which of $\phi$ and $\theta$ denotes the geographical latitude, and whether $\theta=0$ corresponds to the equator, or to the north pole.
To sum it up: There is no general rule. Given the desired "orientation" as well as  a parametrization of the surface $S$ you have to determine "geometrically", i.e. by checking the resulting ${\bf f}_{.1}\times{\bf f}_{.2}$ at a convenient point of the parameter domain whether you have picked the "right orientation".
A: Example 1
Use a parameterisation to find the flux of $F = (3xy, 0, -z)$ outward (normal away from the z-axis) through the cone $z = 9\sqrt{x^2 + y^2}$, $0 \le z \le 9$.

We will use cylindrical coordinates to parameterise the cone.
The formula for the flux through a surface is $\iint_S F(S(\rho, \theta)) \cdot \overrightarrow{N} dt$ where $S(\rho, \theta)$ is the parameterised surface and $\overrightarrow{N} =  \dfrac{\partial{S}}{\partial{\rho}} \times \dfrac{\partial{S}}{\partial{\theta}}$.
The parameterised surface is $S(\rho, \theta) = (\rho \cos(\theta), \rho \sin(\theta), 9\rho) \forall \rho \in [0, 1], \theta \in [0, 2\pi]$.
$\dfrac{\partial{S}}{\partial{\rho}} = (\cos(\theta), \sin(\theta), 9)$
$\dfrac{\partial{S}}{\partial{\theta}} = (-\rho \sin(\theta), \rho \cos(\theta), 0)$
$\overrightarrow{N} = \dfrac{\partial{S}}{\partial{\rho}} \times \dfrac{\partial{S}}{\partial{\theta}} = (-9\rho \cos(\theta), -9\rho \sin(\theta), \rho)$
Therefore, the other normal vector is $- \overrightarrow{N} = (9\rho \cos(\theta), 9\rho \sin(\theta), -\rho)$
$F(S(\rho, \theta)) = (3\rho^2 \cos(\theta)\sin(\theta), 0, -9\rho)$
For any vector representing only radius, the outward-pointing normal vector to such a vector will be in the same direction (assuming, of course, the context of a closed surface).
A vector representing only radius is $S(1, 0) = (1, 0, 9)$
The normal vectors to this point on the surface are
$\overrightarrow{N}(1, 0) = (-9, 0, 1)$ and $-\overrightarrow{N}(1, 0) = (9, 0, -1)$.
Therefore, the correct (outward-pointing) normal vector is $- \overrightarrow{N} = (9\rho \cos(\theta), 9\rho \sin(\theta), -\rho)$, since this is the vector with the same direction as a vector representing only the radius of the surface at the same point $(1, 0)$
$\int_{\theta = 0}^{2\pi} \int_{\rho = 0}^{1} (3\rho^2 \cos(\theta) \sin(\theta), 0, -9\rho) \cdot (9\rho \cos(\theta), 9\rho \sin(\theta), -\rho) d\rho d\theta = 6\pi$

Example 2
Use a parameterisation to find the flux upward across the portion of the plane $x + y + z = 6a$ that lies above the square $0 \le x \le a$, $0 \le y \le a$ in the xy-plane.

To parameterise this surface in terms of $u$ and $v$, let $u = x$, $y = v$, and find a representation in terms of $z$.
$\therefore z = 6a - u - v \ \forall \ 0 \le u \le a$, $0 \le v \le a$
$\therefore S(u, v) = (u, v, 6a - u - v) \ \forall \ 0 \le u \le a$, $0 \le v \le a$
$\dfrac{\partial{S}}{\partial{u}} = (1, 0, -1)$
$\dfrac{\partial{S}}{\partial{v}} = (0, 1, -1)$
$\dfrac{\partial{S}}{\partial{u}} \times \dfrac{\partial{S}}{\partial{v}} = (1, 1, 1)$
$\overrightarrow{S}(1, 0) = (1, 0, 6a - 1)$
$(1, 1, 1)$ has the same direction as $\overrightarrow{S}(1, 0)$.
Therefore, $(1, 1, 1)$ is the outward-pointing normal vector.
$F(S(u, v)) = (4uv, 4v(6a - u - v), 4u(6a - u - v))$
Flux $= \int_{u = 0}^{a} \int_{v = 0}^{a} (4uv, 4v(6a - u - v), 4u(6a - u - v)) \cdot (1, 1, 1) dv du = \dfrac{61a^4}{3}$
A: $\newcommand{\bbx}[1]{\,\bbox[8px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$

Use $\ds{x}$ and $\ds{y}$ as a parameters. Namely, $\ds{\vec{r} = x\,\hat{x} + y\,\hat{y} + 9\root{x^{2} + y^{2}}\hat{z}}$

Then,
\begin{align}
\partiald{\vec{r}}{x}\times\partiald{\vec{r}}{y} & =
\pars{\hat{x}\ +\  {9x \over \root{x^{2} + y^{2}}}\hat{z}}\times
\pars{\hat{y}\ +\  {9y \over \root{x^{2} + y^{2}}}\hat{z}}
\\[5mm] & =
-\,{9x \over \root{x^{2} + y^{2}}}\,\hat{x} -
{9y \over \root{x^{2} + y^{2}}}\,\hat{y} + \hat{z}
\end{align}

\begin{align}
\int\vec{F}\pars{\vec{r}}\cdot\dd\vec{\mrm{S}} & =
\iint\pars{3xy\,\hat{x} - {9 \over \root{x^{2} + y^{2}}}\hat{z}}\cdot
\pars{{9x \over \root{x^{2} + y^{2}}}\,\hat{x} +
{9y \over \root{x^{2} + y^{2}}}\,\hat{y} - \hat{z}}\dd x\,\dd y
\\[5mm] & =
\iint_{\root{x^{2} + y^{2}} < 1}{27x^{2}y + 9 \over \root{x^{2} + y^{2}}}
\,\dd x\,\dd y =
9\int_{0}^{2\pi}\int_{0}^{1}{1\over r}\,r\,\dd r\,\dd\phi = \bbx{\ds{18\pi}}
\end{align}
