Divergence, gradient and differentiation - radial irrotational fluid flow Given a fluid with the steady spherically symmetric flow with only radial velocity $\vec v(r)$. 
We need to evaluate $ \vec v \cdot \nabla \vec v $. 
From vector calculus $$ \vec v \cdot \nabla( \vec v) = \frac{1}{2} grad(v^2) -\vec v \times rot(\vec v)$$
As fluid is irrotational the second term vanishes, so
$$ \vec v \cdot \nabla( \vec v) = \frac{1}{2} grad(v^2)$$
The expression on the right hand side is $$\frac{1}{2} grad(v^2)=\vec v \cdot \frac {\partial \vec v}{\partial r}$$ because we evaluate gradient by simply differentiating $v^2(r)$ with respect to $r$. But term on the left hand side can be written using divirgence in spherical coordinates as
$$ \vec v \cdot div( \vec v) = \vec v \cdot \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 v)$$ (where $v=v_r$ is radial projection of the velocity, scalar value) which after differentiating the product $r^2 v(r)$ gives
$$ \vec v \cdot div( \vec v) =  \vec v \cdot (\frac{2 v}{r}+\frac{\partial v}{\partial r}) = \frac{2\vec v \cdot v}{r} +\vec v \cdot \frac{\partial v}{\partial r}$$
So, there is additional term $\frac{2\vec v \cdot v}{r} $ comparing to the above for $\frac{1}{2}grad(v^2)$ which contradicts to the first expression. Where did I miss something? Many thanks!
 A: OPs fundamental issue is I believe the misunderstanding of the symbol $\vec{\nabla}\vec{v}$ (note: there is no dot product) as a divergence $\text{div}(\vec{v})$. The former is a rank two tensor the latter a scalar.
Going line by line, we have the first equation is a vector equation (vector dotted into tensor equals gradient of scalar):
$$\vec{v}\cdot (\vec{\nabla} \vec{v}) = \frac{1}{2}\vec{\nabla}(v^2)$$
the RHS of this can be written explicitly knowing that $v=v(r)$ only and using the product rule and the fact that for spherical symmetry $\vec{\nabla} = \frac{\partial}{\partial r}\hat{r}$ to give
$$\vec{v}\cdot (\vec{\nabla} \vec{v}) = v \frac{\partial v}{\partial r} \hat{r}$$
note that again, both sides are vectors. The $v$s appearing on the RHS are the scalar magnitude, the direction is in the unit vector $\hat{r}$. There is no dot product on the right.
OPs question is ultimately, why the formula for divergence in spherical coordinates:
$$\vec{\nabla} \cdot \vec{A} = \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 A_r)+ \dots$$
cannot be applied in this case and the short answer is that there are no divergences being taken. There might be a question as to why this nearly works, but I think there's only so many derivative like combinations of a single component vector that this might just be chance.

Remark: I've used $\vec{\nabla}$ (nabla/del) notation throughout this answer since then you can work out the tensor rank of any equation by counting the number of arrows and subtracting two for each dot product present. The index notation also makes this clear since $[\vec{v}\cdot (\vec{\nabla} \vec{v})]_i = v_j\partial_jv_i \neq v_i \partial_j v_j$.
A: Misunderstanding comes from notation. Actually
$$ \vec v \cdot \nabla \vec v \ne \vec v \cdot div(\vec v)$$
The equations in original question just lead to the identity
$$ \vec v \cdot \nabla \vec v = \vec v \cdot div(\vec v)- \frac{2 v^2}{r}\hat{r}$$
for this particular case. The question is closed.
A: I guess I will be the third person to confirm that $$\vec{v} \cdot \nabla \vec{v} \neq \nabla \cdot \vec{v} \, \vec{v} $$ In another type of notation, you can see that
$$\vec{v} \cdot \nabla \vec{v} = \vec{v}\cdot \text{grad}\, \vec{v}$$ while 
$$\nabla \cdot \vec{v} \, \vec{v} = \text{div}(\vec{v})\, \vec{v}$$
According to me, it is more accurate to write this as 
$$\big(\vec{v} \cdot \nabla\big) \, \vec{v}$$
This is actually the term that appears usually in the kinetic part of the various fluid dynamics equations, e.g. Euler equations
\begin{align}
\frac{\partial \vec{v}}{\partial t} + \big(\vec{v} \cdot \nabla\big) \, \vec{v}  &= - \frac{1}{\rho_0}\, \nabla p + \vec{g}\\
\nabla \cdot \vec{v} &= 0
\end{align} 
If your vector field is as you describe it, then I strongly recommend switching to spherical coordinates, because then the vector field $\vec{v}$ simplifies drastically.  Then you would simply work with one component only. Indeed, if your vector field is spherically symmetric and points only in the radial direction, then there exists a function $f(|\vec{r}|) = f(r)$ such that
$$\vec{v} = f(r)\, \vec{r}$$ where $\vec{r} = \begin{bmatrix}x\\y\\z\end{bmatrix}$  and $r = \sqrt{x^2 + y^2 + z^2}$. So basically, in cartesian coordinates 
$$\vec{v} \cdot \nabla = f(r) \left(x \frac{\partial}{\partial x} + y \frac{\partial}{\partial y} + z \frac{\partial}{\partial z}\right)$$ Therefore, in spherical coordinates $(r, \theta, \phi)$, there exists a function $g(r)$ such that 
 $$\vec{v} \cdot \nabla = g(r) \frac{\partial}{\partial r}$$ To calculate it, simply recall that $(\nabla r) \cdot \vec{v} = (\nabla r) \cdot \left(g(r) \frac{\partial}{\partial r}\right) = g(r) (\nabla r) \cdot \left(\frac{\partial}{\partial r}\right) = g(r)$, so when you calculate $(\nabla r) \cdot \vec{v}$ in cartesian coordinates, you get 
\begin{align}
g(r) = (\nabla r) \cdot \vec{v} &= (\nabla r) \cdot \big(f(r)\, \vec{r}\big) = f(r) \Big((\nabla r) \cdot \vec{r}\Big)\\
&=f(r) \, \left(  \Big(\nabla \, \sqrt{x^2 + y^2 + z^2}\Big) \cdot \Big(x \frac{\partial}{\partial x} + y \frac{\partial}{\partial y} + z \frac{\partial}{\partial z}\Big)\right)\\
&= f(r) \, \frac{1}{2 \sqrt{x^2 + y^2 + z^2}}\, 2 \,[x, y, z] \begin{bmatrix}x\\y\\z\end{bmatrix} \\
&= \frac{f(r)}{r}\, \big(x^2 + y^2 + z^2\big) = \frac{f(r)}{r}\, r^2\\
& = r\,f(r)
\end{align} Therefore, the vector field $\vec{v}$ ins polar coordinates is $\vec{v} =\begin{bmatrix}r\, f(r)\\0\\0\end{bmatrix}$ and so
$$\vec{v} \cdot \nabla = r\, f(r) \, \frac{\partial}{\partial r}$$ so 
$$(\vec{v} \cdot \nabla) \, \vec{v} = r\, f(r) \, \frac{\partial}{\partial r} \begin{bmatrix}r\, f(r)\\0\\0\end{bmatrix} = \begin{bmatrix}\Big( r\,f(r)^2 + r^2 \, f(r)\, f'(r) \Big) \\0\\0\end{bmatrix}  $$ 
Observe that from the calculation before $x \frac{\partial}{\partial x} + y \frac{\partial}{\partial y} + z \frac{\partial}{\partial z} = \vec{r} \cdot \nabla = r\, \frac{\partial}{\partial r}$
so going back to cartesian coordinates
$$(\vec{v} \cdot \nabla) \, \vec{v} = \Big( f(r)^2 + r \, f(r)\, f'(r) \Big) \begin{bmatrix} x \\y\\z\end{bmatrix} =  \Big( f(r)^2 + r \, f(r)\, f'(r) \Big)\, \vec{r} $$ where as before $\vec{r} = \begin{bmatrix} x \\y\\z\end{bmatrix}$ and $r = \sqrt{x^2 + y^2 + z^2}$.
