How should one present curl and divergence in an undergraduate multivariable calculus class?

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.]

As far as explaining the formulas for div and curl, you should be able to do this starting with the definitions given in the Wikpedia articles by taking the corresponding integrals on rectangles and boxes. These definitions have fairly clear physical meanings, at least if your students are comfortable with line and surface integrals.

As far as explaining d^2 = 0: curl(grad f) = 0 because the line integral of a gradient about small circles is zero, so gradients can't curl. (In other words, if a vector field has nonzero curl at some p, you wouldn't be able to define a consistent potential about some small closed contour around p.) And div(curl F) = 0 because the surface integral of a curl about small spheres is zero, so curls can't diverge (that is, flow). These interpretations get used all the time in applications of Stokes' theorem to physics.

This is perhaps a crude (and certainly non-rigorous) explanation, but it's always how I've thought of motivating it.

Let $F = (F_1, F_2, F_3)$ denote a vector field in $\mathbb{R}^3$, and write $\text{curl}\ F = (G_1, G_2, G_3)$. We would like a situation where $G_1$ describes the "instantaneous" rotation of $F$ about the $x$-axis, $G_2$ the rotation about the $y$-axis, and $G_3$ the rotation about the $z$-axis.

So let's think of vector fields which do just that. Three simple (linear!) ones which come to mind are $$H_1(x,y,z) = (0, -z, y)$$ $$H_2(x,y,z) = (z, 0, -x)$$ $$H_3(x,y,z) = (-y, x, 0)$$ So in order to measure how much $F$ rotates about, say, the $z$-axis, it makes sense to look at something that compares how similar $F$ is to $H_3$. The dot product $F(x,y,z) \cdot H_3(x,y,z)$ seems reasonable, which is precisely $-yF_1(x,y,z) + xF_2(x,y,z).$

This suggests that defining $$G_1(x,y,z) \approx -zF_2(x,y,z) + yF_3(x,y,z)$$ $$G_2(x,y,z) \approx zF_1(x,y,z) - xF_3(x,y,z)$$ $$G_3(x,y,z) \approx -yF_1(x,y,z) + xF_2(x,y,z)$$ might give something close to what we want. But this is a very crude way to measure "instantaneous" rotation -- in fact, one might say it's a sort of linear approximation. Thus, we are led to replacing the linear terms with their corresponding derivations: $$G_1(x,y,z) = -\frac{\partial}{\partial z}F_2 + \frac{\partial}{\partial y}F_3$$ $$G_2(x,y,z) = \frac{\partial}{\partial z}F_1 - \frac{\partial}{\partial x}F_3$$ $$G_3(x,y,z) = -\frac{\partial}{\partial y}F_1 + \frac{\partial}{\partial x}F_2,$$ which is precisely the curl.

This heuristic also works with divergence, but instead consider $(H_1, H_2, H_3) = (x,y,z)$.