Since your question is how to gain an intuitive understanding of what the two equations $\nabla \cdot \nabla \mathbf{F} = 0$ and $\nabla \times \nabla f=\mathbf{0}$ mean, I will discuss the equations in a simple setting. The simpler, more intuitive setting I am thinking of is a discrete setting: the honeycomb lattice, shown below.
$\hskip 2 in$
Let me explain how vector calculus works on this lattice. A scalar function takes on values at the discrete vertices of the lattice.
A vector field takes on values on the edges of the lattice. This makes sense because vectors should represent movement from one point to another.
The gradient of a function $f$ has, on each edge, a magnitude equal to the difference between $f$ on the two ends of the edge, and the direction of the gradient, of course, is from the lower value to the higher value.
Now when you take the curl you get a function defined on each face (that is, hexagon) of the lattice. This makes sense because the curl is supposed to tell you how the vector field is going around in a "circle". In this case, the "circle" is really a hexagon. Let's see an example of taking the curl of a vector field.
In the above figure, the values of the vector field $\mathbf{F}$ are shown in red. Here I have taken the counter-clockwise direction to be positive, so a negative value means the vector actually points in the clockwise direction. The curl can be found by adding the values as you move counter-clockwise along the hexagon. So the value of the curl at the hexagon shown in the figure is $4$.
Now lets see why the curl of the gradient has to be zero.
Above, I have shown the values of the function in black at each vertex. The gradient is shown in red at each edge as in the previous picture. We find in this case that the curl is zero. The reason is simple. When you are taking the curl in this case, you are adding up the differences in the function values as you move around the hexagon. But since you end up where you started, the sum of the changes must be zero.
Now let's think about why the divergence of the curl must be zero. To take the divergence, we must move to three dimensions. So let's think about the soccer ball shape shown below.
As before, functions are defined on vertices, vectors are defined on edges, the curl is defined on faces. The new thing is that the divergence is defined on the volume enclosed by the soccer ball.
To see why the divergence of the curl is zero, lets just look at what is going on in the two hexagons closest to us.
In this case, we have a vector field $\mathbf{F}$, which can take on arbitrary values on the edges. I have picked random values for the sake of illustration. Again, the curl defined on each hexagon is the sum of $\mathbf{F}$ as you go around the hexagon counterclockwise. The counter-clockwise direction is important here. Notice the edge where $\mathbf{F}=4$ (the one shown in a box) is included in the calculation of the curl at both hexagons. In the lower hexagon, it has a positive sign, since this $\mathbf{F}$ does point in the counterclockwise direction as you go around the lower hexagon. However, it has a negative sign in the upper hexagon, since there this $\mathbf{F}$ points in the clockwise direction. Now since the divergence is the sum the values over all the faces, and the edge's contribution to the lower face is the opposite of the contribution to the upper face, we see that this edge has zero contribution to the divergence. But every single edge on the soccer ball lies on two faces, so each edge gives zero contribution to the divergence, and thus the divergence must be zero.
Now I should probably say a few words about how this relates to the situation of normal euclidean space. The situation is really very analogous. You can measure the curl of a vector field by taking its line integral around small circles. However, in the case of a gradient, the line integral tells you the total change in the function as you go around the circle. Since you end where you begin, the total change must be zero. Since these integrals must all be zero for the gradient, the curl of a gradient must be zero.
The divergence can be measured by integrating the field that goes through a small sphere. If you have a non-zero vector on the surface, then it will tend to create an outward pointing curl on its left, but an inward pointing curl on its right. These two contributions to the divergence cancel out, so you wind up with zero divergence.
The object $\langle \vec{u}, \nabla\rangle$ is an operator defined by $\langle \vec{u}, \nabla\rangle = u_x \frac{\partial}{\partial x} + u_y \frac{\partial}{\partial y} + u_z \frac{\partial}{\partial z}$. It acts on whatever is on the right. In your case it acts on the vector field $\vec{u}(x,y,z,t)$. The action is on each component of $\vec{u}$ separately, so $\langle \vec{u}, \nabla\rangle \vec{u}$ is a vector field:
$$\langle \vec{u}, \nabla\rangle \vec{u} = \left( {\begin{array}{*{20}{c}}
\langle \vec{u}, \nabla\rangle u_x\\
\langle \vec{u}, \nabla\rangle u_y \\
\langle \vec{u}, \nabla\rangle u_z
\end{array}} \right)
=\left( {\begin{array}{*{20}{c}}
u_x \frac{\partial u_x}{\partial x} + u_y \frac{\partial u_x}{\partial y} + u_z \frac{\partial u_x}{\partial z}\\
u_x \frac{\partial u_y}{\partial x} + u_y \frac{\partial u_y}{\partial y} + u_z \frac{\partial u_y}{\partial z} \\
u_x \frac{\partial u_z}{\partial x} + u_y \frac{\partial u_z}{\partial y} + u_z \frac{\partial u_z}{\partial z}
\end{array}}\right).$$
I think the easiest way of understanding where this operator comes from is by imagining following a small parcel of material as it moves along with the flow. For instance, how does the temperature of the parcel change over the short time $dt$?
The parcel is located at position $\vec{r}_1=\left(\matrix{x\\y\\z}\right)$ at time $t_1$ and after a time $dt$ has traveled to $\vec{r}_2$. Since $\vec{u}$ is the parcel's velocity, at time $t_2=t_1+dt$ the parcel is at $\vec{r}_2=\vec{r}_1+\vec{u}dt$.
The temperature field $T(\vec{r},t)$ is a scalar function of $x,y,z$ and $t$. The difference in temperature of the parcel between $t_1$ and $t_2$ is (to lowest order in $dt$),
$$\begin{align}
T(\vec{r}_2, t_2) - T(\vec{r}_1, t_1) &= T(\vec{r}_1+\vec{u}dt, t_1+dt) - T(\vec{r}_1, t_1)\\
&=T(x+u_xdt,y+u_ydt,z+u_zdt,t_1+dt)-T(x,y,z,t_1)\\
&=T(x,y,z,t) + \frac{\partial T}{\partial x} u_x dt + \frac{\partial T}{\partial y} u_y dt + \frac{\partial T}{\partial z} u_z dt + \frac{\partial T}{\partial t} dt - T(x,y,z,t)\\
&= \left( \frac{\partial T}{\partial x} u_x + \frac{\partial T}{\partial y} u_y + \frac{\partial T}{\partial z} u_z + \frac{\partial T}{\partial t} \right) dt\\
&= \left[ \left(\langle \vec{u}, \nabla \rangle + \frac{\partial}{\partial t} \right)T \right] dt.
\end{align}
$$
The rate of change of the temperature of the parcel is the above divided by $dt$. If you wanted to get the rate of change of the $x$-component of the parcel's velocity you would replace $T$ by $u_x$, and similarly for the $y$ and $z$ components.
Best Answer
Your second line of thinking is closer to the truth.
It is correct that the gradient of any harmonic function (i.e., a function satisfying $\nabla^2 \phi = 0$) yields a divergence-free and curl-free vector field. Moreover, the potential for any uniform field $\vec{E}$ is $\phi = - \vec{E} \cdot \vec{r}$, which (it can be shown) is a harmonic function. So it is true that all uniform vector fields can be derived from a harmonic potential.
However, there are also non-uniform vector fields that derive from harmonic potentials. For example, $\phi = \frac12 k(x^2 - y^2)$ is a harmonic function, and it leads to a vector field of $\vec{E} = k(-x \hat{\imath} + y \hat{\jmath})$, which is not uniform.
The reason this seems to be at odds with the Helmholtz Theorem is that the Helmholtz Theorem requires a set of boundary conditions to uniquely determine a vector field, along with its divergence and curl. In particular, if we require that $\vec{E} \to 0$ as $|\vec{r}| \to \infty$, then the Helmholtz theorem tells us that the only divergence- and curl-free vector field that satisfies this is $\vec{E} = 0$ everywhere.1 As we have seen above, however, there are many other divergence- and curl-free vector fields which are possible if we do not require that $\vec{E} \to 0$ at infinity.
1 Note that Laplace's equation also requires a set of boundary conditions to uniquely determine a solution. In particular, if we require that $\nabla^2 \phi = 0$ and $\phi \to \text{constant}$ as $|\vec{r}| \to \infty$, then we must have $\phi = \text{constant}$ everywhere. This is consistent with the fact that the only divergence- and curl-free vector field that vanishes at infinity is the zero vector field.