As far as I know, there is only a really good elementary proof for shapes like cubes, spheres, or polyhedra. Proving this for general regions requires a bit more work.
The simplest way to prove the statement you are asking in the general case is to first prove Gauss's Divergence Theorem. Then for any volume $V$ (with piecewise-smooth boundary), we have
$$
\int\!\!\!\!\int_{S(V)} (\textbf{F}\cdot\textbf{n})\,dS \;=\; \int\!\!\!\!\int\!\!\!\!\int_V (\mathrm{div}\;\textbf{F})\,dV
$$
As $V$ shrinks toward the point $p$, the integral on the right approaches $\mathrm{div}\;\textbf{F}(p)$ multiplied by the volume of $V$. We can make this precise using the following lemma:
Lemma. Let $f\colon\mathbb{R}^n\to\mathbb{R}$ be a continuous function, and let $p\in\mathbb{R}^n$. Then for every $\epsilon>0$, there exists a $\delta>0$ with the following property: if $V$ is a region in $\mathbb{R}^n$ with positive volume and $V$ lies entirely within a distance $\delta$ of the point $p$, then
$$
\left|f(p) - \frac{1}{\mathrm{vol}(V)}\int\!\!\!\!\int\!\!\!\!\int_V f\,dV\right| < \epsilon.
$$
Proof: Let $\epsilon>0$. Since $f$ is continuous, there exists a $\delta > 0$ so that
$$
\|x-p\|<\delta \qquad\Rightarrow\qquad |f(x)-f(p)|<\epsilon
$$
for all $x\in\mathbb{R}^n$. Let $V$ be a region in $\mathbb{R}^n$, and suppose that $V$ lies entirely within $\delta$ of the point $p$. Then $f(p) - \epsilon < f(x) < f(p)+\epsilon$ for all $x\in V$, so
$$
[f(p) - \epsilon]\,\mathrm{vol}(V) \;\leq\; \int\!\!\!\!\int\!\!\!\!\int_V f\,dV \;\leq\; [f(p) + \epsilon]\,\mathrm{vol}(V)
$$
and the result follows.$\qquad\square$
Using this theorem and Gauss's Divergence Theorem, we immediately get:
Theorem. Let $\textbf{F}\colon\mathbb{R}^n\to\mathbb{R}^n$ be a continuously differentiable vector field, and let $p\in\mathbb{R}^n$. Then for every $\epsilon>0$, there exists a $\delta>0$ with the following property: if $V$ is a region in $\mathbb{R}^n$ with positive volume and a piecewise-smooth boundary and $V$ lies entirely within a distance $\delta$ of the point $p$, then
$$
\left|\mathrm{div}\;\textbf{F}(p) - \frac{1}{\mathrm{vol}(V)}\int\!\!\!\!\int_{S(V)} (\textbf{F}\cdot\textbf{n})\,dS\right| < \epsilon.\tag*{$\square$}
$$
It is possible to restate this theorem using sequences and limits. For the following theorem, if $V$ is a region in $\mathbb{R}^n$, define
$$
\mathrm{rad}(V,p) \;=\; \sup\{\|x-p\| \mid x\in V\}.
$$
That is, $\mathrm{rad}(V,p)$ is the radius of the smallest sphere centered at $p$ that contains $V$. Then
Theorem. Let $\textbf{F}\colon\mathbb{R}^n\to\mathbb{R}^n$ be a continuously differentiable vector field, and let $p\in\mathbb{R}^n$. Let $V_1,V_2,V_3,\ldots$ be a sequence of regions in $\mathbb{R}^n$, all with positive volume and piecewise-smooth boundary, and suppose that
$$
\lim_{n\to\infty} \mathrm{rad}(V,p) \;=\; 0.
$$
Then
$$
\lim_{n\to\infty} \frac{1}{\mathrm{vol}(V_n)}\int\!\!\!\!\int_{S(V_n)} (\textbf{F}\cdot\textbf{n})\,dS \;=\; \mathrm{div}\;\textbf{F}(p).\tag*{$\square$}
$$
