Here's a sketch of the relation between div-grad-curl and the de Rham complex, in case you might find it useful.
The first thing to realise is that the div-grad-curl story is inextricably linked to calculus in a three-dimensional euclidean space. This is not surprising if you consider that this stuff used to go by the name of "vector calculus" at a time when a physicist's definition of a vector was "a quantity with both magnitude and direction". Hence the inner product is essential part of the baggage as is the three-dimensionality (in the guise of the cross product of vectors).
In three-dimensional euclidean space you have the inner product and the cross product and this allows you to write the de Rham sequence in terms of div, grad and curl as follows:
$$ \matrix{ \Omega^0 & \stackrel{d}{\longrightarrow} & \Omega^1 & \stackrel{d}{\longrightarrow} & \Omega^2 & \stackrel{d}{\longrightarrow} & \Omega^3 \cr
\uparrow & & \uparrow & & \uparrow & & \uparrow \cr
\Omega^0 & \stackrel{\mathrm{grad}}{\longrightarrow} & \mathcal{X} & \stackrel{\mathrm{curl}}{\longrightarrow} & \mathcal{X} & \stackrel{\mathrm{div}}{\longrightarrow} & \Omega^0 \cr}$$
where $\mathcal{X}$ stands for vector fields and the vertical maps are, from left to right, the following isomorphisms:
- the identity: $f \mapsto f$
- the musical isomorphism $X \mapsto \langle X, -\rangle$
- $X \mapsto \omega$, where $\omega(Y,Z) = \langle X, Y \times Z \rangle$
- $f \mapsto f \mathrm{dvol}$, where $\mathrm{dvol}(X,Y,Z) = \langle X, Y \times Z\rangle$
up to perhaps a sign here and there that I'm too lazy to chase.
The beauty of this is that, first of all, the two vector calculus identities $\mathrm{div} \circ \mathrm{curl} = 0$ and $\mathrm{curl} \circ \mathrm{grad} = 0$ are now subsumed simply in $d^2 = 0$, and that whereas div, grad, curl are trapped in three-dimensional euclidean space, the de Rham complex exists in any differentiable manifold without any extra structure. We teach the language of differential forms to our undergraduates in Edinburgh in their third year and this is one way to motivate it.
As for the integral theorems, I always found Spivak's Calculus on manifolds to be a pretty good book.
Another answer mentioned Gravitation by Misner, Thorne and Wheeler. Personally I found their treatment of differential forms very confusing when I was a student. I'm happier with the idea of a dual vector space than I am with the "milk crates" they draw to illustrate differential forms. Wald's book on General Relativity had, to my mind, a much nicer treatment of this subject.
To me, the explanation for the appearance of div, grad and curl in physical equations is in their invariance properties.
Physics undergrads are taught (aren't they?) Galileo's principle that physical laws should be invariant under inertial coordinate changes. So take a first-order differential operator $D$, mapping 3-vector fields to 3-vector fields. If it's to appear in any general physical equation, it must commute with with translations (and therefore have constant coefficients) and also with rotations. Just by considering rotations about the 3 coordinate axes, you can then check that $D$ is a multiple of curl.
If I want to devise a "physical" operator which has the same invariance property - and therefore equals curl, up to a factor - I'd try something like "the mean angular velocity of particles uniformly distributed on a very small sphere centred at $\mathbf{x}$, as they are carried along by the vector field." (This is manifestly invariant, but not manifestly a differential operator!)
[Here I should admit that, having occasionally tried, I've never convinced more than a fraction of a calculus class that it's possible to understand something in terms of the properties it satisfies rather than in terms of a formula. That's unsurprising, perhaps: it's not an obvious idea, and it's entirely absent from the standard textbooks.]
Best Answer
Let ${\bf K}$ be a vector field in the neighbourhood of ${\bf p}\in{\mathbb R}^n$, and let ${\bf X}$ and ${\bf Y}$ be two tangent vectors at ${\bf p}$. These two vectors span a parallelogram $P$ with one vertex at ${\bf p}$. The "circulation" of ${\bf K}$ around $P$ computes to $$ \int_{\partial P}{\bf K}\cdot \mathrm{d}{\bf x}= (L\,{\bf X})\cdot{\bf Y}- (L\,{\bf Y})\cdot{\bf X} + o(|P|^2) $$ with $L:=\mathrm{d}{\bf K}({\bf p})$ and $|P|:= \mathrm{diam}(P)$. It follows that there is a certain skew bilinear function ${\rm Rot}{\bf K}({\bf p}):T_{\bf p}\times T_{\bf p}\to{\mathbb R}$ with $$ \int_{\partial P}{\bf K}\cdot \mathrm{d}{\bf x}={\rm Rot}{\bf K}({\bf p})({\bf X},{\bf Y})+ o(|P|^2) \quad (|P|\to 0). $$ In the case $n=3$ the bilinear form ${\rm Rot}$ can be represented by the vector ${\rm curl}{\bf K}$ in the form $$ {\rm Rot}{\bf K}({\bf p})({\bf X},{\bf Y}) = {\rm curl}{\bf K}({\bf p})\cdot({\bf X}\times{\bf Y}). $$