$$\rho\frac{D \vec V}{Dt}=-\nabla p+ \mu \nabla^{2} \vec V+\rho g$$
In the Navier-Stokes equations there's this term $ \mu \nabla^{2} \vec V $. I don't really understand what this means. What is the physical meaning of the Laplacian of the velocity vector field? I've seen some explanations of this, but they seem complicated. Could you just give an intuitive explanation of what it means?
Since you are asking for some intuition, forget about the other terms in the Navier-Stokes apart from:
$$ \partial_t \bf v = \mu \nabla^2v $$
All the other terms of the full equation will do their job (convection, external forces...), but this simplification highlights the pure effect of the Laplacian. Its effect is not modified by the presence of all the other terms, just all the processes driven by each term will happen at the same time in the full equation.
The above equation is the simplest version of the diffusion equation, not different from the heat equation.
The name already gives you the idea of what this Laplacian term is doing: it will smooth out the gradients of the velocity field and will tend to make the velocity field more and more uniform as time passes by: this is the diffusion of momentum in a viscous fluid (the heat equation is the diffusion of internal kinetic energy at the molecular level).
Assume a steady and homogeneous flow in the whole space (this is a trivial solution of the above equation) and perturb it in some limited region of space: after a while, the solution has relaxed back to the original, unperturbed configuration. This is the typical behaviour of a dissipative (i.e., "viscous") fluid (let the water in your glass relax to a stationary state, perturb it, and observe it to relax back to the initial zero-velocity state: the mechanical energy of the perturbed flow at the macroscopic scale is now in some internal degrees of freedom of the liquid).
This is the physical effect of the Laplacian term, but why does the Laplacian implement the "diffusion" process? My intuition (that applies to this case as well as to the usual heat equation) is this:
More explicitly:
$$ \partial_t \text{(valley)} = \nabla^2 \text{(valley)} >0 \qquad \text{(i.e., the valleys go up)} $$
$$ \partial_t \text{(mount top)} = \nabla^2 \text{(mount top)} <0 \qquad \text{(i.e., the peaks go down)} $$
In short, it's an internal friction of the fluid. The friction is proportional to how fast the fluid velocity changes sideways to the velocity. For a piece of fluid following the flow, $$ \mathbf{F}_{\text{friction}} = \mu \iint_{\text{surface of piece of fluid}} (\mathbf{n}\cdot\nabla)\mathbf{v} \, dS, $$ where $\mathbf{n}$ is the unit normal. (Note that $\mathbf{n}\cdot\mathbf{v}=0$ on the surface of the piece of fluid; no fluid enters or leaves the piece through the surface.)
By a generalization of the divergence theorem, the friction force can be written$$ \mathbf{F}_{\text{friction}} = \mu \iiint_{\text{internal of piece of fluid}} \nabla^2\mathbf{v} \, dV, $$ which is why $\mu \nabla^2\mathbf{v}$ occurs in the Navier-Stokes equation.
I'll derive the Navier-Stokes equation, but focusing on the viscosity term, that is the critical one. Indeed I'll write only a justification for the equation: the price to pay to offer a simple road to the Navier-Stokes equation is that it will not be a proof, but only a reasonable argument. I hope that I (or someone else) will improve these steps in future, keeping the road to $\eta \nabla^2 \bar{v}$ simple, but more satisfying. Before reading this answer it could be very useful (if not essential) reading the answer I wrote on the Cauchy equation and on the Euler equation. Note: I will use no bar for scalars, one bar for vectors (hat for versors) and two bars for tensors (for brevity, when I will say "tensor" I will always mean "rank 2 tensor", but of course scalars and vectors too are tensors).
Navier-Stokes is similar to the Euler equation but it describes a fluid with viscosity (indeed we find the Euler equation taking $\eta \to 0$). In evaluating forces in the interior of the fluid we take into account not only forces due to pressure, but also of forces that are dependent on how fluid moves. It is intuitive that these "extra" forces are zero if fluid moves with rigid motion (think to Newton bucket), and are bigger if motion try very hard to destroy rigid motion (think to a turbine), but transforming this vague intuition into equations is very hard and in what follow I'll do some heuristic steps (I found this by myself almost alone because Cengel-Cimbala doesn't report proofs and Granger is often too hard and obscure: I hope that writing this "work done" will prevent to others my headaches). It would be better if these heuristics steps were linked with some microscopic interpretation or some direct experimental evidence (in this sense, many headache problems are left open here: are there simple road to fix these steps?).
Navier-Stokes equation is (be not afraid, I'll calmly explain each term) \begin{equation} \rho \frac{D \bar{v}}{D t} = \rho \bar{f} - \nabla p + \eta \nabla^2 \bar{v} \tag{1} \end{equation} where $\rho$ is density (constant), $\bar{v}$ velocity field, $p$ is pressure, $\eta$ is viscosity (see below) and $\bar{f}$ is the vector described in my S.E. answer about the Cauchy equation. The operator $\frac{D}{Dt}$ is the material derivative $\frac{\partial }{\partial t} + \bar{v} \cdot \nabla \otimes$ (see my S.E. answer about Cauchy equation). Finally, $\nabla^2 \bar{v}$ is by definition the vector whose components are laplacian of components: $(\nabla^2 v_x , \nabla^2 v_y , \nabla^2 v_z)$.
Observe that $\bar{v}(\bar{r},t) = v_x(\bar{r},t)\hat{x} + v_y(\bar{r},t)\hat{y} + v_z(\bar{r},t)\hat{z}$ so omitting many dependencies (when not written they are $(\bar{r},t)$) we can write \begin{equation} \bar{v}(\bar{r}+d\bar{r},t) = \left( v_x + \frac{\partial v_x}{\partial x}dx + \frac{\partial v_x}{\partial y}dy + \frac{\partial v_x}{\partial z}dz \right) \hat{x} + \left( v_y + \frac{\partial v_y}{\partial x}dx + \frac{\partial v_y}{\partial y}dy + \frac{\partial v_y}{\partial z}dz \right) \hat{y} + \left( v_z + \frac{\partial v_z}{\partial x}dx + \frac{\partial v_z}{\partial y}dy + \frac{\partial v_z}{\partial z}dz \right) \hat{z} \end{equation} And so \begin{equation} \bar{v}(\bar{r}+d\bar{r},t) = \bar{v} (\bar{r},t) + \left( \frac{\partial v_x}{\partial x}dx + \frac{\partial v_x}{\partial y}dy + \frac{\partial v_x}{\partial z}dz \right) \hat{x} + \left( \frac{\partial v_y}{\partial x}dx + \frac{\partial v_y}{\partial y}dy + \frac{\partial v_y}{\partial z}dz \right) \hat{y} + \left( \frac{\partial v_z}{\partial x}dx + \frac{\partial v_z}{\partial y}dy + \frac{\partial v_z}{\partial z}dz \right) \hat{z} \end{equation} Exploiting the notation introduced in my S.E. answer about Cauchy equation, for example with \begin{equation} \nabla \otimes \bar{v} = \begin{pmatrix} \frac{\partial v_x}{\partial x} & \frac{\partial v_y}{\partial x} & \frac{\partial v_z}{\partial x} \\ \frac{\partial v_x}{\partial y} & \frac{\partial v_y}{\partial y} & \frac{\partial v_z}{\partial y} \\ \frac{\partial v_x}{\partial z} & \frac{\partial v_y}{\partial z} & \frac{\partial v_z}{\partial z} \end{pmatrix} \end{equation} We could also write $\bar{v}(\bar{r}+d\bar{r},t) = \bar{v} (\bar{r},t) + \nabla \otimes \bar{v} \cdot d \bar{r}$, a nicer equation but it is not essential do that now.
We have found the different speed of two nearby point (at a given time $t$), now we can check that the different speed we have found, can be written in this way too (Granger call this "Helmholtz theorem") \begin{equation} \bar{v}(\bar{r}+d\bar{r},t) = \bar{v} (\bar{r},t) + \bar{\omega} \times d\bar{r} + \bar{\bar{\epsilon}} \cdot d \bar{r} \tag{2} \end{equation} where \begin{equation} \bar{\omega} = \frac{1}{2} \nabla \times \bar{v} \end{equation} and \begin{equation} \bar{\bar{\epsilon}} = \begin{pmatrix} \frac{\partial v_x}{\partial x} & \frac{1}{2} \left( \frac{\partial v_x}{\partial y} + \frac{\partial v_y}{\partial x} \right) & \frac{1}{2} \left( \frac{\partial v_x}{\partial z} + \frac{\partial v_z}{\partial x} \right) \\ \frac{1}{2} \left( \frac{\partial v_x}{\partial y} + \frac{\partial v_y}{\partial x} \right) & \frac{\partial v_y}{\partial y} & \frac{1}{2} \left( \frac{\partial v_y}{\partial z} + \frac{\partial v_z}{\partial y} \right) \\ \frac{1}{2} \left( \frac{\partial v_x}{\partial z} + \frac{\partial v_z}{\partial x} \right) & \frac{1}{2} \left( \frac{\partial v_y}{\partial z} + \frac{\partial v_z}{\partial y} \right) & \frac{\partial v_z}{\partial z} \end{pmatrix} \tag{3} \end{equation} The check is long but not difficult, in principle.
This symmetric tensor can be written in a nicer way as $\bar{\bar{\epsilon}} = \frac{\nabla \otimes \bar{v} + (\nabla \otimes \bar{v})^T}{2}$, as you can check. In literature you may find $\bar{\bar{\epsilon}}$ here or there written in that way (but $\otimes$ is omitted) and it is called strain rate tensor, but these shorter notations are not important now.
The tensor $\bar{\bar{\epsilon}}$ is the mathematical term that describes the speed of "deformation" of the fluid, I mean that this term describes how much the fluid deviates from rigid motion. Why do I claim that? Because the fundamental equation of rigid motion says that each time $\exists ! \,\, \bar{\omega} : \bar{v}_A = \bar{v}_B + \bar{\omega} \times (\bar{r}_A - \bar{r}_B) \, \, \forall \, A,B$ (as consequence of the constraint $|\bar{r}_A - \bar{r}_B|$=constant $\forall \, A,B$), so looking at (2) I see that the claim is meaningful. Of course a fluid doesn't have a relaxed state to which it tends to return (while a solid has) and speaking about deformation of a fluid can sound wrong, but we are speaking about speed of deformation, and this is meaningful. We have that the difference of speed between $\bar{r}$ and $\bar{r}+d\bar{r}$ (at a given time $t$) has a part that is linked to how much fluid deviates from rigid motion, and this part is (infinitesimal because it is proportional to distances between the two points, who are close by hipotesis) \begin{equation} d\tilde{v} = \bar{\bar{\epsilon}} \cdot d \bar{r} \end{equation} Please note that $d\tilde{v}$ is an infinitesimal vector, function of $(\bar{r},t,d\bar{r})$. It is linked to the speed of deformation but we don't see how such a vector can be used in determining viscous forces: we have to work on it. To evaluate viscous forces we would like to have a vector that is bigger when speed of deformation is bigger (and zero in case of rigid motion), but that is finite and that is function of $(\bar{r},t,\hat{n})$, i.e. of position, time and direction. It is reasonable that a such vector multiplied by $dA$ would give the force exerted on the atoms lying near one side of the imaginary surface $d\bar{A}=\hat{n}dA$ inside the fluid, by the atoms near the other side. A such vector can be found by dividing by $dr$: \begin{equation} \frac{d\tilde{v}}{dr} = \bar{\bar{\epsilon}} \cdot \hat{n} \end{equation} In conclusion, our hypotesis is that viscous forces acting on the "lower" side of a surface $d\bar{A}$ inside the fluid are proportional to $dA$ and to $\bar{\bar{\epsilon}} \cdot \hat{n}$ with $\hat{n}=\frac{d\bar{A}}{dA}=\frac{d\bar{r}}{dr}$ (I move in direction $d\bar{A}$ by a lenght $dr$, and the I normalize finding $\hat{n}$, $\bar{\bar{\epsilon}} \cdot \hat{n} dr$ is the "not rigid" contribution to the different speed between starting point $\bar{r}$ and end point $\bar{r}+d\bar{r}$). In short, we assume that for some constant $k$, the force exerted on the "lower" side of $d\bar{A}$ due to the presence of viscosity is \begin{equation} d \bar{F} = k \bar{\bar{\epsilon}} \cdot d\bar{A} \tag{4} \end{equation} i.e. (look to definition of stress tensor in my answer on Cauchy equation) we assume that the stress tensor that describes viscous forces is \begin{equation} \bar{\bar{\sigma}}_{visc} = k \bar{\bar{\epsilon}} \tag{5} \end{equation}
Now we will show how the $k$ constant introduced is linked to viscosity. Viscosity is defined in a bit artificial way: with two infinite parallel plates sliding (you will find things about that everywhere in books and web). You can't create a such configuration in laboratory, so it looks at first glance too abstract, but as pointed out by Cengel-Cimbala there are practical flows for which this "Couette flow" idealization is a very good approximation, this allows experimental confirmation of newtonian law of viscosity. For example using a long coaxial cylinders viscometer ($R_i$ and $R_o$ inner and outer radius, length $L$, torque applied $N$ and angular speed $\omega$) we find experimentally that for a given fluid the quantity $\frac{N(R_o-R_i)}{2\pi\omega R_i^3 L}$ is always the same constant $\eta$, and it is the same constant that we find in Poiseuille law $\frac{dm}{dt} = \frac{\rho \pi R^4}{8 \eta} \frac{dp}{dz}$. This shows that the newtonian definition of viscosity works. So the abstract, artificial way in which viscosity is defined is acceptable.
Let's so focus on mental experiment that define viscosity: let's think to two plates that slide on each other (separated by the fluid), and suppose that the surface $d\bar{A}$ is in the lower plate oriented upwards, as in figure
Axes are at rest with the lower plate. $V\hat{y}$ is the speed of the upper plates. $\hat{n}=\hat{z}$ and exploiting (4) and (3) we see that \begin{equation} d\bar{F} = k \left( \frac{1}{2} \left( \frac{\partial v_x}{\partial z}+\frac{\partial v_z}{\partial x} \right) , \frac{1}{2} \left( \frac{\partial v_y}{\partial z}+\frac{\partial v_z}{\partial y}\right) , \frac{\partial v_z}{\partial z} \right) dA = \left( 0 , \frac{k}{2} \frac{\partial v_y}{\partial z} , 0\right) dA \end{equation} because velocity field has only $y$ component and varies only in $z$ direction. Effectively we see by simmetry that velocity field must be $\bar{v}=\left( 0,\frac{V}{h} z , 0 \right)$, so the stress acting on $d\bar{A}$ is (don't confuse stress tensor $\bar{\bar{\sigma}}$ with the stress vector $\bar{\sigma}=\frac{d\bar{F}}{dA}$, in our case orizontal) $\bar{\sigma} = \frac{kV}{2h}\hat{y}$. But by definition of viscosity its magnitude is ${\sigma} = \eta \frac{dv}{dz}$, so we have $\bar{\sigma} = \eta \frac{V}{h} \hat{y}$, and we conclude that $k=2\eta$. Exploiting (5) we get the stress tensor that describes viscous forces: \begin{equation} \bar{\bar{\sigma}}_{visc} = 2 \eta \bar{\bar{\epsilon}} \tag{6} \end{equation} Note that this tensor is $\bar{\bar{0}}$ if $\eta \to 0$ and/or if $\bar{v}$ is uniform (not necessarily constant in time, think to the zero viscous forces of water in a vessel on the ground of an elevator).
We know (see answer about Euler equation) that for inviscid fluid (zero viscosity) the stress tensor is $\bar{\bar{\sigma}}=-p\bar{\bar{I}}$, and we know that the stress tensor that describes viscous forces is given by (6). The simpler solution is to assume that taking account of both pressure and viscosity, the stress tensor is simply the sum of the two: \begin{equation} \bar{\bar{\sigma}}=-p\bar{\bar{I}} + 2 \eta \bar{\bar{\epsilon}} \tag{7} \end{equation} Often in studying nature we start with simple hypotheses (think to Hooke law, Fourier law of heat, Planck hypothesis, etc.) and then we test them to check if they work and we try to go in deep in understanding the meaning. Of course I suppose this try of writing stress tensor can be better justified, but like Cengel-Cimbala I'll simply claim claim to generalize the $-p\bar{\bar{I}}$ adding another stress tensor that take into account viscous forces (but unlike what they do, I added some arguments to try to justify this term that describes the deviation from rigid motion). Even if I don't find it satisfactory, I find it reasonable.
Please note that $-p\bar{\bar{I}}$ gives only normal stresses (see my S.E. answer about Euler equation) while $2 \eta \bar{\bar{\epsilon}} $ gives both normal and tangential stresses.
In the final of my S.E. answer about Cauchy equation I wrote it in this way \begin{equation} \frac{D \bar{v}}{D t}=\frac{1}{\rho} \nabla \cdot \bar{\bar{\sigma}} + \bar{f} \tag{8} \end{equation} If stress tensor is writable as in (7), we have then \begin{equation} \rho \frac{D \bar{v}}{D t}= \left( \frac{\partial}{\partial x} \quad \frac{\partial}{\partial y} \quad \frac{\partial}{\partial z} \right) \begin{pmatrix} -p & 0 & 0 \\ 0 & -p & 0 \\ 0 & 0 & -p \end{pmatrix} + 2 \eta \left( \frac{\partial}{\partial x} \quad \frac{\partial}{\partial y} \quad \frac{\partial}{\partial z} \right) \begin{pmatrix} \epsilon_{xx} & \epsilon_{xy} & \epsilon_{xz} \\ \epsilon_{xy} & \epsilon_{yy} & \epsilon_{yz} \\ \epsilon_{xz} & \epsilon_{yz} & \epsilon_{zz} \end{pmatrix} + \rho \bar{f} \end{equation} Let's focus on the $x$ component, the first big term is clearly $-\frac{\partial p}{\partial x}$. The other can be calculated by exploiting (3): \begin{equation} 2 \eta \left( \frac{\partial \epsilon_{xx}}{\partial x} + \frac{\partial \epsilon_{xy}}{\partial y} + \frac{\partial \epsilon_{xz}}{\partial z} \right) = \eta \left( 2 \frac{\partial^2 v_x}{ \partial x^2} + \frac{\partial^2 v_x}{ \partial y^2} + \frac{\partial^2 v_y}{ \partial x \partial y} + \frac{\partial^2 v_x}{ \partial z^2} + \frac{\partial^2 v_z}{ \partial x \partial z} \right) \end{equation} The term in parenthesis can be written as $\frac{\partial^2 v_x}{ \partial x^2} + \frac{\partial^2 v_x}{ \partial y^2}+ \frac{\partial^2 v_x}{ \partial z^2} + \frac{\partial}{\partial x} \left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} + \frac{\partial v_z}{\partial z} \right) $, i.e. as $\nabla^2 v_x + \frac{\partial}{\partial x} (\nabla \cdot \bar{v}) = \nabla^2 v_x$ because the fluid is incompressible so $\nabla \cdot \bar{v}=0$. We conclude that $x$ component of the equation (8) must be \begin{equation} \rho \left[ \frac{D \bar{v}}{D t} \right]_x = - \frac{\partial p}{\partial x} + \eta \nabla^2 v_x + \rho f_x \end{equation} The other components will be similar, so after having defined $\nabla^2 \bar{v}$ as we did after (1), we got Navier-Stokes equation, which is nothing but the Cauchy equation for momentum conservation after having chosen (7) as stress tensor (please note that if we consider $\nabla \cdot \bar{v} = 0$ due the continuity equation, the number of unknown $v_x$, $v_y$, $v_z$ and $p$ is equal to the number of independent equations, but it is not surprising that will not be a picnic exploiting these four equations to explore problems).
