\(\)
Divergence measures how much the flow is expanding. It singles out the stretching component of motion.
AKA as the GaussGreen theorem. The 3D analogue to Green’s theorem for flux.
Definition
If \(S\) is a closed surface, enclosing a region \(D\), oriented with \(\hat n\) outwards and \(\vec F\) defined and differentiable everywhere in \(D\), then $$ \newcommand{oiint}{\subset\!\!\supset\kern1.65em\iint} \shaded{ \oiint_S\vec F\cdot d\vec S = \iiint_D\rm{div}(\vec F)\,dV } \label{eq:divthm} $$ The circle in the double integral reminds us that it must be a closed surface. The normal vectors always point outwards.
Let \(\vec F\) be a vector field, where \(P,Q,R\) are functions of \(x,y,z\). $$ \newcommand{pdv}[2]{\frac{\partial #1}{\partial #2}} \shaded{ \begin{align*} \rm{div}\left(\vec F\right) &= \rm{div}\left(P\hat\imath+Q\hat\jmath+R\hat k\right) \\ &= \pdv{P}{x}+\pdv{Q}{y}+\pdv{R}{z} \\ &= P_x + Q_y + R_z \end{align*} } \nonumber $$
The divergence measures how much the flow is expanding. If you take a region of space, the total amount of water that flows out of it is the total amount of sources that you have in there minus the sinks.
\(\nabla\)notation
\(\nabla\) “del” is a symbolic notation for the operator $$ \newcommand{pdv}[2]{\frac{\partial #1}{\partial #2}} \nabla = \left\langle\pdv{}{x}, \pdv{}{y}, \pdv{}{z}\right\rangle \nonumber $$
That makes $$ \newcommand{pdv}[2]{\frac{\partial #1}{\partial #2}} \begin{align*} \nabla f &= \left\langle\pdv{f}{x}, \pdv{f}{y}, \pdv{f}{z}\right\rangle & \text{gradient} \\[0.7em] \nabla\cdot\vec F &= \left\langle\pdv{}{x}, \pdv{}{y}, \pdv{}{z}\right\rangle \cdot \left\langle P,Q,R\right\rangle \\ &= \pdv{P}{x} + \pdv{Q}{y} + \pdv{R}{z} & \text{divergence} \end{align*} $$
This will come in very useful when we talk about curl in space.
Physical interpretation
Let \(\vec F\) be a vector field, where \(P,Q,R\) are functions of \(x,y,z\). $$ \newcommand{oiint}{\subset\!\!\supset\kern1.65em\iint} \newcommand{pdv}[2]{\frac{\partial #1}{\partial #2}} \begin{align*} \oiint_S\vec F\cdot d\vec S &= \iiint_D\rm{div}(\vec F)\,dV \\[0.5em] \Rightarrow \oiint_S\left\langle P,Q,R\right\rangle &= \iiint_D\left(P_x+Q_y+R_z\right)\,dV \end{align*} $$
Claim that $$ \rm{div}\left(\vec F\right) = \text{“source rate”} \nonumber $$
The amount of flux generated per unit volume. Think about an incompressible fluid, such as water. Incompressible fluid flow with velocity \(\vec F\): (given mass occupies a fixed volume),
Proof of divergence theorem
Simplifications
 The dotproduct can be seen a a sum of the different dimensions. We will proof one component: $$ \newcommand{oiint}{\subset\!\!\supset\kern1.65em\iint} \shaded{ \oiint_S\left\langle 0,0,R\right\rangle \cdot \hat n\,dS = \iiint_D R_z\,dV } \label{eq:simple1} $$ To proof the general case, we will prove the same thing for a vector field that has only an \(x\) or \(y\)component, and then sum the identities. $$ \newcommand{oiint}{\subset\!\!\supset\kern1.65em\iint} \begin{align*} \oiint_S\left\langle P,Q,R\right\rangle\hat n\,dS &= \oiint_S\left\langle P,0,0\right\rangle \cdot \hat n\,dS \\ &+ \oiint_S\left\langle 0,Q,0\right\rangle \cdot \hat n\,dS \\ &+ \oiint_S\left\langle 0,0,R\right\rangle \cdot \hat n\,dS \end{align*} $$
 We will start with a vertically simple region. If region \(D\) is vertically simple It lives above the shadow on the \(xy\)plane, and is between two graphs \(z_1(x,y)\) and \(z_2(x,y)\).
$$ \newcommand{oiint}{\subset\!\!\supset\kern1.65em\iint} \oiint_S\left\langle 0,0,R\right\rangle \cdot \hat n\,dS = \iiint_D R_z\,dz \nonumber $$
We will compute the right and lefthand sides of equation \(\eqref{eq:simple1}\), and then compare them.
Righthand side
In righthand side of equation \(\eqref{eq:simple1}\), the range for \(z\) is from the bottom face to the top face. $$ \begin{align*} \iiint_D R_z\,dV &= \iint_U \underline{\int_{z_1(x,y)}^{z_2(x,y)} R_z\,dz}\,dx\,dy \end{align*} $$
When you integrate \(\newcommand{pdv}[2]{\frac{\partial #1}{\partial #2}} R_z=\pdv{R}{z}\) in respect to \(z\), you just get \(R\) back $$ \begin{align} \iiint_D R_z\,dV &= \iint_U \Big[R(x,y,z)\Big]_{z_1(x,y)}^{z=z_2(x,y)}\,dx\,dy \nonumber \\ &= \iint_U \Big[R(x,y,z_2(x,y))R(x,y,z_1(x,y))\Big]\,dx\,dy \label{eq:divthmproofleft} \\ \end{align} $$
Lefthand side
The surface \(S\) of the vertically simple region, consists of the top, bottom and the sides.
We will calculate the flux on each side, and then sum them together $$ \newcommand{oiint}{\subset\!\!\supset\kern1.65em\iint} \begin{align} \oiint_S\left\langle 0,0,R\right\rangle \cdot \hat n\,dS &= \iint_\rm{top} + \iint_\rm{bottom} + \iint_\rm{sides} \label{eq:combine} \end{align} $$
Top
Recall equation \(\eqref{eq:fluxgraph}\) for the flux of a graph
$$ \hat n\,dS =\pm\left\langle f_x,f_y,1\right\rangle\,dx\,dy \nonumber $$
Apply this to the graph of \(z=z_2(x,y)\) $$ \newcommand{pdv}[2]{\frac{\partial #1}{\partial #2}} \hat n\,dS = \left\langle \pdv{z_2}{x}, \pdv{z_2}{y}, 1 \right\rangle\,dx\,dy \nonumber $$
Do the dotproduct with the vector field that only has a \(z\)component and enter it in the double integral \(\eqref{eq:simple1}\) $$ \begin{align*} \iint_\text{top}\left\langle 0,0,R\right\rangle\cdot\hat n\,dS &= \iint_\text{top}\left\langle 0,0,R\right\rangle\cdot \left\langle \pdv{z_2}{x}, \pdv{z_2}{y}, 1 \right\rangle\\ &= \iint_\text{top} R\,dx\,dy \end{align*} $$
\(R\) is \(R(x,y,z)\). That means we can substitute the equation using the graph equation \(z=z_2(x,y)\). As far as the region is concerned, it is above shadow \(U\) $$ \newcommand{pdv}[2]{\frac{\partial #1}{\partial #2}} \begin{align} \iint_\rm{top} \left\langle 0,0,R\right\rangle\cdot\hat n\,dS &= \iint_U R(x,y,z_2(x,y))\,dx\,dy \label{eq:divthmproofright1} \\ \end{align} $$
Bottom
Similarly, apply equation \(\eqref{eq:fluxgraph}\) to the graph of \(z=z_1(x,y)\) $$ \newcommand{pdv}[2]{\frac{\partial #1}{\partial #2}} \hat n\,dS = \left\langle \pdv{z_1}{x}, \pdv{z_1}{y}, 1 \right\rangle \nonumber $$
Be care with the orientation of the normal vector. At the top surface, we want \(\hat n\) pointing up, but at the bottom surface, we want \(\hat n\) pointing down.
Switch the orientation, so the normal vector points downwards $$ \newcommand{oiint}{\subset\!\!\supset\kern1.65em\iint} \newcommand{pdv}[2]{\frac{\partial #1}{\partial #2}} \begin{align*} d\vec S &= \left\langle \pdv{z_1}{x},\pdv{z_1}{y}, 1\right\rangle\,dx\,dy \end{align*} $$
Do the dotproduct with the vector field that only has a \(z\)component and enter it in the double integral \(\eqref{eq:simple1}\) $$ \begin{align*} \iint_\text{bottom}\left\langle 0,0,R\right\rangle\cdot\hat n\,dS &= \iint_\text{bottom}\left\langle 0,0,R\right\rangle\cdot \left\langle\pdv{z_2}{x},\pdv{z_2}{y},1 \right\rangle\\ &= \iint_\text{bottom} R\,dx\,dy \end{align*} $$
\(R\) is \(R(x,y,z)\). That means we can substitute the equation using the graph equation \(z=z_2(x,y)\). As far as the region is concerned, it is above shadow \(U\) $$ \newcommand{pdv}[2]{\frac{\partial #1}{\partial #2}} \begin{align} \iint_\rm{bottom} \left\langle 0,0,R\right\rangle\cdot\hat n\,dS &= \iint_U R(x,y,z_2(x,y))\,dx\,dy \label{eq:divthmproofright2} \\ \end{align} $$
Sides
\(\left\langle 0,0,R\right\rangle\) is tangent to the sides \(\Longrightarrow\) the flux through the sides \(=0\). (that is why we took a vector field with only a \(z\)component.)
In other words: the vector field \(\left\langle 0,0,R\right\rangle\) is parallel to the \(z\)axis, so at the sides there nothing going on.
Together
Combining equation \(\eqref{eq:divthmproofright1}\) and \(\eqref{eq:divthmproofright2}\) using \(\eqref{eq:combine}\) $$ \newcommand{oiint}{\subset\!\!\supset\kern1.65em\iint} \begin{align} \oiint_S\left\langle 0,0,R\right\rangle \cdot \hat n\,dS = \iint_U &R(x,y,z_2(x,y))\,dx\,dy \nonumber \\ + \iint_U &R(x,y,z_1(x,y))\,dx\,dy \label{eq:together} \end{align} $$
Compare left and righthand side
Comparing \(\eqref{eq:divthmproofleft}\) with \(\eqref{eq:together}\) you indeed get the same formula $$ \newcommand{pdv}[2]{\frac{\partial #1}{\partial #2}} \newcommand{oiint}{\subset\!\!\supset\kern1.65em\iint} \left\{ \begin{align*} \oiint_S\left\langle 0,0,R\right\rangle \cdot \hat n\,dS = \iint_U &\Big[R(x,y,z_2(x,y)) \\ &R(x,y,z_1(x,y))\Big]\,dx\,dy \\ \iiint_D R_z\,dV = \iint_U &R(x,y,z_2(x,y))\,dx\,dy \\ + \iint_U &R(x,y,z_1(x,y))\,dx\,dy \end{align*} \right. $$
Therefore, for a vertically simple region $$ \newcommand{oiint}{\subset\!\!\supset\kern1.65em\iint} \shaded{ \oiint_S\left\langle 0,0,R\right\rangle \cdot \hat n\,dS = \iiint_D R_z\,dV } \nonumber $$
Generalization
Lift the simplifications

If \(D\) is not vertically simple, we cut it into vertically simple pieces.
For instance, a solid doughnut is not vertically simple because above and below the hole makes two regions.
However, if we slice it into four pieces, then each piece has a well defined top and bottom side (=vertically simple).
 We could also do the same proof for the \(x\) and \(y\)component. Then sum them together to get the general case.
Example
Redo the earlier example: Find the flux of \(\vec F=\left\langle x,y,z\right\rangle\) through sphere of radius \(a\) centered at the origin. The normal vector of the sphere points out radially from the origin.
Apply the divergence theorem \(\eqref{eq:divthm}\) $$ \newcommand{oiint}{\subset\!\!\supset\kern1.65em\iint} \begin{align*} \oiint_S\vec z\hat k\cdot d\vec S &= \iiint_D\rm{div}(z\hat k)\,dV \\ &= \iiint_D (0+0+1) \,dV \\ &= \underbrace{\iiint_D dV}_{\text{Volume}(D)} = \frac{4}{3}\pi a^3\\ \end{align*} $$