What is an intuitive explanation for $\operatorname{div} \operatorname{curl} F = 0$? I am aware of an intuitive explanation  for $\operatorname{curl} \operatorname{grad} F = 0$ (a block placed on a mountainous frictionless surface will slide to lower ground without spinning), and was wondering if there were a similar explanation for $\operatorname{div} \operatorname{curl} F = 0$.
 A: The curl of a vector field $F$ at a point on an infinitesimal surface $S$ can be thought of as how the vector field "spins" $S$. The divergence, on the other hand, describes how much the vector field "passes through" $S$. If you take the a vector field which spins a surface, how much does it pass through the surface? None.
Of course, this is a highly non-technical explanation.
A: Think of the surface of water, and a two dimensional vector field $\mathbf{F}(x,y)$ which gives its velocity everywhere. Think that you put a needle that floats somewhere on its surface. If the needle is rotating, then it is in a region of non-zer curl, and if the needle (EDIT: a bunch of them) appears to be translating inwards or outwards towards or away from some point, then it is a region of non zero divergence. Finite divergence just means there's like a hose pipe somewhere which is sucking in or spewing out fluid radially.  This is how you can think of the divergence and the curl individually.
For your question, think of a mixed region with non zero divergence and non-zero curl. The direction of the curl would be along the axis of rotation and hence perpendicular to the surface. 
$$\mathbf{\nabla} \times \mathbf{F}(x,y) = g(x,y) \hat{z}$$
$$\nabla \cdot \mathbf{\nabla} \times \mathbf{F}(x,y) =  \frac{\partial g(x,y) }{\partial z} = 0$$
Or intuitively, the curl is perpendicular to the divergence, which is an abuse of notation because divergence is really an operator, so this is as far as the analogy goes.
A: Update: I felt that in my original answer, the intuition got lost in all of the formalism, so I've rewritten it to make it clearer. The core idea is still much the same as Qiaochu's answer on the linked MO question.
First off, I don't think you're going to get an explanation for $\nabla\cdot (\nabla\times \mathbf F) = 0$ that's on quite the same level as your sliding block example, because while the gravitational force is well-known to be the gradient of the negative potential, there aren't any tangible vector fields that are the curl of something. The curl, like the cross product, is a little "unnatural" as it depends on the choice of handedness; classical physics, on the other hand, is independent of handedness, so you can never observe a curl directly, only hidden under an integral or cross product to make the handedness goes away.
Nevertheless, the most mathematically natural argument that $\nabla\cdot (\nabla\times \mathbf F) = 0$ is quite easy to grasp. For a warm-up, let's look at another explanation for why $\nabla\times (\nabla f) = 0$. The curl of any vector field $\mathbf F$ describes how much $\mathbf F$ "spins around", right? So take any closed loop in space; the net curl inside the loop is nothing but how much $\mathbf F$ "circulates" around the loop, i.e. $\oint \mathbf F \cdot d\mathbf x$. But if $\mathbf F$ is a gradient of something, say $f$, then $\mathbf F \cdot d\mathbf x$ is actually how much $f$ is going up or down as you walk along $d\mathbf x$. So if you walk along a closed loop and come back to where you started, the net change in $f$ has to be zero! Or as wzzx's comment states: you can't walk from home to school and back and have gone uphill both ways.
Now we can do the same thing for $\nabla\cdot(\nabla\times \mathbf F)$. Intuitively, the divergence of a vector field $\mathbf G$ measures how much $\mathbf G$ is "spreading out" or "pulling in". In other words, pick any region of space; what does the total divergence of $\mathbf G$ inside it tell you? It tells you exactly how much $\mathbf G$ is flowing out of the surface of the region, i.e. $\oint \mathbf G \cdot d\mathbf A$. But if $\mathbf G$ is the curl of another vector field $\mathbf F$, then $\int \mathbf G \cdot d\mathbf A$ on a surface just measures how much $\mathbf F$ circulates around the boundary of that surface. What's the boundary of a closed surface? Imagine taking a portion of the entire surface and growing it to cover the whole. Its boundary eventually gets smaller and smaller, and then disappears. So there is nothing for $\mathbf F$ to circulate around, and the circulation must be zero.
All I've shown here is that the integrals of $\nabla\times(\nabla f)$ and $\nabla\cdot(\nabla\times\mathbf F)$ are zero on any arbitrary region, but that should be enough to see that (at least for smooth $f$ and $\mathbf F$) their values must be zero everywhere.
A: Formal Intuition
"Formal intuition" is a bit of a contradiction, but look at 
$$\mathbf{A} \cdot (\mathbf{A} \times \mathbf{B})$$
Cross products are perpendicular to the things being crossed, so $\mathbf{A \times B}$ is perpendicular to $\mathbf{A}$. The dot product of perpendicular vectors is zero, so
$$\mathbf{A} \cdot (\mathbf{A} \times \mathbf{B}) = 0$$
Replacing $\mathbf{A}$ with $\nabla$ would give the identity, but I'm not sure if this counts as intuitive.  Note: This isn't supposed to be a proof.  It's only supposed to suggest the identity.
Geometric Intuition
Here, I'll use Gauss's theorem and Stoke's theorem, which themselves have intuitive explanations.
Make a coordinate system like this:

Consider a cubical box.  Make it centered at the origin with sides of length $1$, so the corners are $(.5,.5,.5)$, $(.5,.5,-.5)$, $(.5,-.5,.5)$, etc.  This is for convenience - the area of each side is $1$ and the volume is $1$, so any time I need to divide or multiply by the area of a side or by the volume, nothing changes.
We are looking at $\nabla \cdot (\nabla \times \mathbf{F})$.  Ignoring the stuff inside the parentheses, we're looking at a divergence of something.  Let's call that thing $\mathbf{C}$ since it's a curl.
Gauss' theorem says that the integral of the divergence over the volume of the box is equal to the net flux coming out.  This means that we can find the average divergence over the entire box by thinking only about the faces.
Consider the top face.  To find the flux coming out of the top face, first find the average of $\mathbf{C}$ over that face.  Then project that vector onto a unit normal vector, $\hat{z}$.  This gives a vector $\mathbf{C}_z$ that either points into or out of the box.  If it points out, it contributes positive flux, and if it points in, it contributes negative flux.

If we can show that all six such vectors (one for each face) contribute zero net flux when combined, we know that the average divergence over our box is zero.  How do we get information about the average of $\mathbf{C}$, projected onto a normal vector, over some face?
This is where the curl comes in.  Stokes' theorem says that the average of $\mathbf{C}_z = (\nabla \times \mathbf{F})\cdot\hat{z}$ over the top is equal to the line integral over the edges of the top face of $\mathbf{F}\cdot \hat{r}$, with $\hat{r}$ a unit tangent vector.  In other words, each of the four edges of the top gives some contribution to the average of $\mathbf{C}_z$ on the top.
Grab one particular edge - the one from $(.5, -.5, .5)$ to $(.5,.5,.5)$ will do.  Suppose we traverse it from left to right, and get a positive number, say $4.2$, out of the line integral of $\mathbf{F} \cdot \hat{r}$ along that edge.

The purple vector is the original vector field $\mathbf{F}$.  It gets projected onto $\hat{r}$ to make $F_r$, the purple segment.  This gets integrated along the edge, and contributes to the black vector $\mathbf{C_z}$, which itself contributes to the flux.
The important thing to note is that this $4.2$ value of the line integral along the edge gets used twice in finding the flux because the edge borders on two different faces.  For the top face, going left-to-right on our edge means circling clockwise as viewed from above.  The right-hand rule tells us that the positive $4.2$ of the line integral along our edge contributes positively to a vector $\mathbf{C_z}$ point out of the top.  Since it points up out of the top, it contributes negative flux.  All told, traversing our edge, thinking of it as part of the top of the cube, contributes $-4.2$ to the total flux.
When we look at the front of the cube, going left-to-right on our edge means going clockwise around the front, as viewed from in front of the cube.  The right hand rule means that a positive clockwise line integral around the front creates a vector $\mathbf{C}_x$ pointing into the cube - positive flux.  So this time our edge contributes $+4.2$ to the total flux.

The edge's two contributions to the flux cancel each other because one point outwards and the other points inwards.  There is nothing special about this edge, so all edges will be canceled by their two contributions to the flux.  The total flux for this cube is zero.  That means the average of the divergence over the cube is zero.  There was nothing special about this cube, so the average divergence over any cube is zero, and the divergence itself must simply be zero everywhere.
Physical Intuition
This theorem about the divergence of the curl is responsible for the conservation of charge.  Begin with Maxwell's equation for the curl of the magnetic field
$$\nabla \times \mathbf{B} = \frac{1}{c}\frac{\partial \mathbf{E}}{\partial t} + \frac{4 \pi}{c} \mathbf{J}$$
Take the divergence of both sides.
$$\nabla \cdot (\nabla \times \mathbf{B}) = \frac{1}{c} \frac{\partial}{\partial t} \nabla \cdot \mathbf{E} + \frac{4\pi}{c}\nabla\cdot\mathbf{J}$$
Now use Maxwell's equation for the divergence of the electric field
$$\nabla \cdot \mathbf{E} = 4\pi \rho$$ 
Insert this into the above to get
$$\nabla \cdot (\nabla \times \mathbf{B}) = \frac{4 \pi}{c} \left(\frac{\partial \rho}{\partial t} + \nabla\cdot\mathbf{J}\right)$$
We expect that electric charge will be conserved.  That means the the flux of charge current coming out of some region must be compensated by a loss of charge density there.
$$\nabla \cdot \mathbf{J} = -\frac{\partial \rho}{\partial t}$$
Plugging this into our previous expression gives
$$\nabla \cdot (\nabla \times \mathbf{B}) = \frac{4 \pi}{c} \left(\frac{\partial \rho}{\partial t} -\frac{\partial \rho}{\partial t}\right) = 0$$
This means that if electromagnetic fields obey Maxwell's equations and if electric charge is conserved, then the divergence of the curl of any magnetic field must be zero.  However, I'm not sure how intuitive this is, and also the magnetic field cannot be chosen arbitrarily because it must be divergenceless.
We could repeat the procedure for Maxwell's equation for the curl of the electric field.  There, we'd get
$$\nabla \cdot (\nabla \times \mathbf{E}) = \frac{\partial}{\partial t}\left(\nabla \cdot \mathbf{B}\right) = 0$$
because the magnetic field is divergenceless.  So the absence of magnetic charge implies that the divergence of the curl of all electric fields is zero.  Electric fields, unlike magnetic fields, can be essentially arbitrary as far as I know, and so viewed in a certain light Maxwell's equations imply the sought vector identity. 
A: Here's an explanation I was given while I was still in high school (but not in school, of course!). I didn't fully understand it at the time, but it might be illuminating if you're willing to accept the generalised Stokes theorem. Essentially, it's the claim that
$$\int_{M} d\omega = \int_{\partial M} \omega$$
where $\omega$ is a differential $(n - 1)$-form, $d\omega$ is its exterior derivative, $M$ is an $n$-dimensional smooth manifold with boundary, and $\partial M$ its boundary. (There are some technical conditions which I'm glossing over.) I should probably explain what these things are though: 


*

*A smooth $n$-manifold with boundary is an $n$-dimensional smooth geometric object where every interior point has a neighbourhood which looks like a piece of our familiar Euclidean $n$-space $\mathbb{R}^n$, and every boundary point has a neighbourhood which looks like a piece of $\mathbb{R}^{n-1} \times [0, \infty)$, which is just $\mathbb{R}^n$ chopped in half. For example, the interval $[0, 1]$ is a smooth 1-manifold with boundary, and its boundary is the set $\{ 0, 1 \}$; and the closed unit disc is a 2-manifold with boundary the unit circle.

*Note that the unit circle is a 1-manifold without boundary; this is an instance of the general fact that the boundary of a smooth $n$-manifold with boundary is a smooth $(n-1)$-manifold without boundary.

*A differential $n$-form is essentially what you integrate when doing an $n$-fold integral. For example, in $\displaystyle \int_{-\infty}^{\infty} \exp(-x^2) \, dx$, you are integrating the differential 1-form $\exp(-x^2) \, dx$ over the 1-manifold $\mathbb{R}$, and in $\displaystyle \int_{0}^{2\pi} \int_0^1 r \, dr \, d\theta$, you are integrating the differential 2-form $r \, dr \wedge d\theta$ over the closed unit disc.

*The exterior derivative of a differential $(n-1)$-form is a differential $n$-form. It's simultaneously a generalisation of the gradient, curl, and divergence operators. (More on this later.)


Now, I suppose I should explain what this has to do with div, grad, and curl. Let's start with a function $f$ of 3 real variables $(x, y, z)$. This can be regarded as a differential 0-form, and so it has an exterior derivative $df$. This is essentially the same thing as $\nabla f$; indeed, we can expand $df$ using the chain rule to get
$$df = \frac{\partial f}{\partial x} \, dx + \frac{\partial f}{\partial y} \, dy + \frac{\partial f}{\partial z} \, dz$$
which, neglecting certain subtleties, is the same as $\nabla f$ expanded out in the standard basis. Suppose instead I give you a differential 1-form $F = a \, dx + b \, dy + c \, dz$. Its exterior derivative is
$$dF = \left( \frac{\partial c}{\partial y} - \frac{\partial b}{\partial z} \right) \, dy \wedge dz + \left( \frac{\partial a}{\partial z} - \frac{\partial c}{\partial x} \right) \, dz \wedge dx + \left( \frac{\partial b}{\partial x} - \frac{\partial a}{\partial y} \right) \, dx \wedge dy$$
which, with a little squinting, is seen to be more-or-less the curl $\nabla \times F$. (This is basically the chain rule again, but you might be wondering why there aren't 9 terms in the exterior derivative and where all the minuses are coming from. This is because for any differential 1-form $\alpha$, $\alpha \wedge \alpha = 0$. Applying this to $(dx + dy) \wedge (dx + dy)$, we see that $dx \wedge dy = -dy \wedge dx$ — that is, the wedge product $\wedge$ of differential 1-forms is anticommutative, just like the cross product of two vectors.)
Finally, suppose I have a differential 2-form $A = a \, dy \wedge dz + b \, dz \wedge dx + c \, dx \wedge dy$. Its exterior derivative is
$$dA = \left( \frac{\partial a}{\partial x} + \frac{\partial b}{\partial y} + \frac{\partial c}{\partial z} \right) \, dx \wedge dy \wedge dz$$
which is quite obviously related to the divergence $\nabla \cdot A$.
The generalised Stokes theorem is a vast generalisation of the fundamental theorem of calculus, the Kelvin—Stokes theorem, and the divergence theorem. Indeed, they are all special cases: again, working in 3-dimensional Euclidean space, we get


*

*When $\omega$ is a 0-form, the fundamental theorem of calculus: $$\displaystyle \int_{P}^{Q} \nabla f \cdot d\vec{x} = f(Q) - f(P)$$

*When $\omega$ is a 1-form, the Kelvin—Stokes theorem: $$\displaystyle \int_S \nabla \times \vec{F} \cdot d\vec{S} = \oint_{\partial S} \vec{F} \cdot d\vec{x}$$

*When $\omega$ is a 2-form, the divergence theorem: $$\displaystyle \int_V \nabla \cdot \vec{A} \, dV = \int_{\partial V} \vec{A} \cdot d\vec{S}$$


Of course, exactly as with div, grad, and curl, just by fiddling with the definitions, you can show that for any differential $n$-form $\omega$, the exterior derivative of the exterior derivative is $d^2 \omega = 0$. But this doesn't explain why it should happen, and I think the best explanation is through the generalised Stokes theorem:
$$\int_M d^2 \omega = \int_{\partial M} d\omega = \int_{\partial^2 M} \omega$$
This shows that the fact that $d^2 \omega = 0$ is intimately linked to the fact that $\partial^2 M = \emptyset$, i.e. that the boundary of a smooth manifold itself doesn't have a boundary. Translating this to the classical notation, we immediately have the identities
$$\nabla \times (\nabla f) = 0$$
$$\nabla \cdot (\nabla \times F) = 0$$
as promised.
A: I apologize first if my explanation below had been proposed before because there are so many answers that I have no time to check every comment. I also wish to thank my student in my Electromagnetic Theory class, Yung-Chi Chien, for helpful discussions this year.
To show that div curf F = 0, the simplest way is to expand curl F and then div curl in the Cartesian coordinate system. Changing the order of partial differentiation will cancel out all of the expanded terms.
Here is a "proof" without using specific coordinate systems. Consider a surface integral of curl F over an arbitrary closed surface, say, a spherical surface S.  Due to the divergence theorem, that surface integral equals the volume integral of div curl F.
At the same time, let the sphere S be cut into two.  The surface integral of curl F over the upper hemisphere equals a line integral of F along a circle determined by the boundary curve C (a circle here) of the two hemispheres.  The line integral on C may be, for example, counterclockwise. Similarly, the surface integral of curl F over the lower hemisphere is again equal to a line integral on C but integrating clockwise.  These two line integrals on C thus cancel each other.  Hence the surface integral over the whole S will be zero, and the volume integral of div curl F is zero, too.  Since we can change the radius and position of the sphere at will, div curl F must be zero.
curl grad phi can be proved similarly, and many people have pointed out clear interpretations.
