I'm trying to follow an explanation of the intuition behind 2-D divergence (article here) using Green's Theorem in the normal form, but I am having trouble understanding the last (and probably simplest) step. Any help would be appreciated.
Green's Theorem gives that the flux on a vector field $\mathbf{\vec F}$ over a closed curve C is equal to the double integral over the enclosed region of C of the divergence of $\mathbf{\vec F}$ (provided the region is continuously differentiable), namely, $\oint_C\mathbf{\vec F} \cdot \mathbf{\hat n}\,ds = \iint_R \frac{\partial M}{\partial x} + \frac{\partial N}{\partial y}\,dA$, where $\mathbf{\vec F} = M(x,y) \mathbf{\hat i} + N(x,y) \mathbf{\hat j}$ and represents a velocity field (fluid flow field).
The intuition proceeds to explain the integrand $\frac{\partial M}{\partial x} + \frac{\partial N}{\partial y}\,dA$ as follows. Take a very small rectangle in the region R, where the bottom left of the rectangle has coordinates $(x,y)$, and which has a width and height of $\Delta x$ and $\Delta y$, respectively. Calculating the flux across the boundary of this rectangle (going in the counterclockwise direction) gives approximately the integrand in question.
To see this, consider the top of the rectangle. The flux across this piece will be approximately $\mathbf{\vec F}(x,y+\Delta y) \cdot \mathbf{\vec j}\,\Delta x = N(x,y+\Delta y)\Delta x$. We've assumed that the rectangle is small enough to where $\mathbf{\vec F}$ is constant across the top of the rectangle, when in reality it would change, so instead of taking the integral at each point to calculate net flux, we just use the dot product to calculate the area of the parallelogram.
Similarly, the flux across the bottom of the rectangle is approximately $\mathbf{\vec F}(x,y) \cdot -\mathbf{\vec j}\,\Delta x = -N(x,y) \Delta x$, again assuming that $\mathbf{\vec F}(x+\Delta x, y)$, and everywhere in between, is the same as $\mathbf{\vec F}(x,y)$.
Adding the two fluxes, we get that the flux is approximately $(N(x,y+\Delta y) - N(x,y))\Delta x$. Since $\frac{N(x,y+\Delta y) - N(x,y)}{\Delta y} \approx \frac{\partial N}{\partial y}$, we can rewrite the flux as approximately $(\frac{\partial N}{\partial y}\Delta y)\Delta x$.
Similar reasoning shows that the flux across the left and right sides of the rectangle is approximately $(\frac{\partial M}{\partial x}\Delta x)\Delta y$. Therefore, the total flux is approximately the sum of these two effects, namely, $(\frac{\partial M}{\partial x} + \frac{\partial N}{\partial y})\Delta x \Delta y$, which gives the integrand.
Because of conservation of matter, fluid that is flowing out of the rectangle must be accounted for by fluid flowing in. Therefore, the flux across this rectangle is equal to the source rate (the rate at which fluid is coming in) of the rectangle, so we can say that $(\frac{\partial M}{\partial x} + \frac{\partial N}{\partial y})\Delta x \Delta y$ approximates the source rate.
Here is the part that I don't understand:
Source rate of rectangle $\approx (\frac{\partial M}{\partial x} + \frac{\partial N}{\partial y})\Delta x \Delta y$
Dividing by the area and "passing to the limit" gives
Source rate at (x,y) $= \frac{\partial M}{\partial x} + \frac{\partial N}{\partial y} = div \,\mathbf{\vec F}$, and therefore, the 2-D divergence represents the source rate at a point $(x,y)$.
I'm not following exactly how this division by the area works, probably because the left sides of the equations are English instead of symbols. Can someone explain mathematically what is going on in this last step?
Bonus question: we're assuming $\mathbf{\vec F}$ is constant across each side of the rectangle in order to derive the total flux out of the rectangle. For instance, $\mathbf{\vec F}$ over the left side of the rectangle is $\mathbf{\vec F}(x,y)$, meaning that even at the top of the left side, $(x,y+\Delta y)$, we're taking $\mathbf{\vec F}$ to be $\mathbf{\vec F}(x,y)$. However, on the top side of the rectangle, we're assuming $\mathbf{\vec F}$ is $\mathbf{\vec F}(x,y+\Delta y)$ everywhere. This seems to contradict what we just assumed for the left side of the rectangle, since we took $\mathbf{\vec F}$ at $(x,y+\Delta y)$ to be $\mathbf{\vec F}(x,y)$. Why can't we just assume the rectangle is small enough so that $\mathbf{\vec F}$ is constant everywhere? Of course, if we do that then we get the total flux is 0. Why do some approximations like this work, and others don't?
Well, you want to consider divergence at the point, so you take the limit as the area goes to $0$. The small rectangle you considered had area $\Delta x \Delta y$. Now, the total divergence over the area is approximately $(\frac{\partial M}{\partial x} + \frac{\partial N}{\partial y}) \Delta x \Delta y$, so to consider it at a point, you divide by the area, leaving you with $\frac{\partial M}{\partial x} + \frac{\partial N}{\partial y}$.
Your bonus question is good. This is an intuitive explanation so it's hard to make it rigorous, but the idea is that $F$ varies over the rectangle so let's consider how it varies over the $x$ and $y$ directions separately, keeping the other variable constant, a technique used often in multivariable calculus. The fact that you choose $x = x$ when considering the flux across the top to be $F(x,y+\Delta y) \cdot j$ is actually not particularly relevant. In fact, I would argue that you should use $F(x + \frac{\Delta x}{2}, y + \Delta y)$ to use an average $x$-value over the rectangle. Regardless, they used intuitive approximations, and it lacks mathematical rigour. As for why some approximations work and others don't, it's because some assumptions are reasonable and some are not. It's reasonable to decouple the changes in the $x$ and $y$ directions while keeping the other variable constant, but it would be less reasonable to assume that the value of $F$ over the rectangle is always constant.