7.1. Fluid Dynamics¶

7.1.1. Stress-Energy Tensor¶

In general, the stress energy tensor is the flux of momentum $p^\mu$ over the surface $x^\nu$ . It is a machine that contains a knowledge of the energy density, momentum density and stress as measured by any observer of the event.

Imagine a (small) box in the spacetime. Then the observer with a 4-velocity $u^\mu$ measures the density of 4-momentum $\d p^\alpha\over\d V$ in his frame as:

${\d p^\alpha\over\d V} = -T^\alpha{}_\beta u^\beta$

and the energy density that he measures is:

$\rho = {E\over V} = -{u^\alpha p_\alpha \over V} = - u^\alpha {\d p_\alpha\over\d V} = u^\alpha T_{\alpha\beta} u^\beta$

One can also obtain the stress energy tensor from the Lagrangian $\L=\L(\eta_\rho, \partial_\nu \eta_\rho, x^\nu)$ by combining the Euler-Lagrange equations

${ \partial \L\over\partial \eta_\rho} - \partial_\nu\left( { \partial \L\over\partial (\partial_\nu \eta_\rho)} \right) =0$

with the total derivative ${\d \L\over \d x^\mu}$ :

${\d \L\over \d x^\mu} = {\partial\L\over\partial\eta_\rho} \partial_\mu \eta_\rho + { \partial \L\over\partial (\partial_\nu \eta_\rho)} \partial_\mu\partial_\nu\eta_\rho + \partial_\mu\L = = \partial_\nu\left( { \partial \L\over\partial (\partial_\nu \eta_\rho)} \right) \partial_\mu \eta_\rho + { \partial \L\over\partial (\partial_\nu \eta_\rho)} \partial_\nu\partial_\mu\eta_\rho + \partial_\mu\L = = \partial_\nu\left( { \partial \L\over\partial (\partial_\nu \eta_\rho)} \partial_\mu \eta_\rho \right) + \partial_\mu\L$

or

$\partial_\nu\left( { \partial \L\over\partial (\partial_\nu \eta_\rho)} \partial_\mu \eta_\rho -\L \delta_\mu{}^\nu \right) + \partial_\mu\L =0$

This can be written as:

$\partial_\nu T_\mu{}^\nu + f_\mu = 0$

where

$T_\mu{}^\nu = { \partial \L\over\partial (\partial_\nu \eta_\rho)} \partial_\mu \eta_\rho -\L \delta_\mu{}^\nu f_\mu = \partial_\mu\L$

The Navier-Stokes equations can be derived from the conservation law:

$\partial_\nu T^{\mu\nu} + f^\mu = 0$

To obtain some Lagrangian (and action) for the perfect fluid, so that we can derive the stress energy tensor $T^{\mu\nu}$ from that, is not trivial, see for example arXiv:gr-qc/9304026. One has to take into account the equation of state and incorporate the particle number conservation $\nabla_\mu(nu^\mu)=0$ and no entropy exchange $\nabla_\mu(nsu^\mu)=0$ constraints.

The equation of continuity follows from the conservation of the baryon number — the volume $V$ that contains certain number of baryons can change, but the total number of baryons $nV$ must remain constant:

${\d (nV)\over\d\tau} = 0 {\d n\over\d\tau}V + n{\d V\over\d\tau} = 0 u^\alpha (\partial_\alpha n)V + n(\partial_\alpha u^\alpha) V = 0 \partial_\alpha (n u^\alpha) = 0$

Perfect Fluids¶

Perfect fluids have no heat conduction ( $T^{i0} = T^{0i} = 0$ ) and no viscosity ( $T^{ij} = p\one$ ), so in the comoving frame:

$T^{\alpha\beta} = \diag(\rho c^2, p, p, p) = \left(\rho+{p\over c^2}\right)u^\alpha u^\beta + p g^{\alpha\beta}$

where in the comoving frame we have $g^{\mu\nu} = \diag(-1, 1, 1, 1)$ , $u^0=c$ and $u^i=0$ , but $\partial_\alpha U^i\neq0$ . $p$ is the pressure with units $[p] =\rm N\,m^{-2}=kg\,m^{-1}\,s^{-2}$ (then $[{p\over c^2}] =\rm kg\,m^{-3}$ ), $\rho$ is the rest mass density with units $[\rho] =\rm kg\,m^{-3}$ , and $\rho c^2$ is the energy density with units $[\rho c^2] =\rm kg\,m^{-1}\,s^{-2}$ .

The last equation is a tensor equation so it holds in any frame. Let’s write the components explicitly:

$T^{00} = \left(\rho+{p\over c^2}\right)u^0u^0 - p = \left(\rho+{p\over c^2}\right)c^2 \gamma^2 - p = \left(\rho c^2+p\left(1-{1\over\gamma^2}\right)\right) \gamma^2 = \left(\rho c^2+p {v^2\over c^2}\right) \gamma^2 T^{0i} = T^{i0} = \left(\rho+{p\over c^2}\right)u^0u^i = \left(\rho+{p\over c^2}\right) c v^i \gamma^2 = {1\over c}\left(\rho c^2+p\right) v^i \gamma^2 T^{ij} = \left(\rho+{p\over c^2}\right) u^iu^j + p \delta^{ij} = \left(\rho+{p\over c^2}\right) v^iv^j\gamma^2 + p \delta^{ij}$

We now use the conservation of the stress energy tensor and the conservation of the number of particles:

(7.1.1.1)¶ $\partial_\nu T^{\mu\nu} = 0$

(7.1.1.2)¶ $\partial_\mu(nu^\mu) = 0$

The equation (7.1.1.2) gives:

$\partial_t (n\gamma) + \partial_i(n v^i \gamma) = 0$

(7.1.1.3)¶ $\partial_t (n m\gamma) + \partial_i(n m v^i \gamma) = 0$

(7.1.1.4)¶ $\partial_t (n m c^2\gamma) + \partial_i(n m c^2 v^i \gamma) = 0$

The equation (7.1.1.1) gives for $\mu=0$ :

$\partial_\nu T^{0\nu} = 0 \partial_0 T^{00} + \partial_i T^{0i} = 0 \partial_t\left({1\over c}\left(\rho c^2 + p {v^2\over c^2}\right) \gamma^2\right) + \partial_i\left({1\over c}\left(\rho c^2 + p\right) v^i \gamma^2\right) = 0$

(7.1.1.5)¶ $\partial_t\left(\left(\rho c^2 + p {v^2\over c^2}\right) \gamma^2\right) + \partial_i\left(\left(\rho c^2 + p\right) v^i \gamma^2\right) = 0$

We now substract the equation (7.1.1.4) from (7.1.1.5):

$\partial_t\left(\left(\rho c^2\gamma - n m c^2 + p {v^2\over c^2} \gamma\right) \gamma\right) + \partial_i\left(\left(\rho c^2\gamma -n m c^2 + p \gamma\right) v^i \gamma\right) = 0$

We define the nonrelativistic energy as:

$E = \rho c^2\gamma -n m c^2 = \half \rho v^2 + (\rho - nm)c^2 + O\left(v^4\over c^2\right)$

so it contains the kinetic plus internal energies. We substitute back into (7.1.1.5):

(7.1.1.6)¶ $\partial_t\left(\left(E + p {v^2\over c^2} \gamma\right) \gamma\right) + \partial_i\left(\left(E + p \gamma\right) v^i \gamma\right) = 0$

This is the relativistic equation for the energy. Substituting $nm = \rho\gamma - {E\over c^2}$ into (7.1.1.3):

(7.1.1.7)¶ $\partial_t\left(\rho\gamma^2 - {E\gamma\over c^2}\right) + \partial_i\left(\left(\rho\gamma^2 - {E\gamma\over c^2}\right) v^i \right) = 0$

The equation (7.1.1.1) for $\mu=i$ gives:

$\partial_\nu T^{i\nu} = 0 \partial_0 T^{i0} + \partial_j T^{ij} = 0 \partial_t \left({1\over c^2}\left(\rho c^2 + p \right) v^i\gamma^2\right) + \partial_j \left( \left(\rho+{p\over c^2}\right) v^iv^j\gamma^2 + p \delta^{ij} \right) = 0$

(7.1.1.8)¶ $\partial_t \left(\left(\rho + {p\over c^2} \right) v^i\gamma^2\right) + \partial_j \left( \left(\rho+{p\over c^2}\right) v^iv^j\gamma^2 + p \delta^{ij} \right) = 0$

This is the momentum equation. The equations (7.1.1.7), (7.1.1.8) and (7.1.1.6) are the correct relativistic equations for the perfect fluid (no approximations were done). We can take either (7.1.1.7) or (7.1.1.5) as the equation of continuity (both give the same nonrelativistic equation of continuity). Their Newtonian limit is obtained by $c\to\infty$ (which implies $\gamma\to1$ ):

$\partial_t \rho + \partial_i(\rho v^i) = 0 \partial_t \left(\rho v^i\right) + \partial_j \left( \rho v^iv^j + p \delta^{ij} \right) = 0 \partial_t E + \partial_j\left(v^j\left(E + p \right)\right) = 0$

those are the Euler equations, also sometimes written as:

(7.1.1.9)¶ ${\partial \rho\over \partial t} + \nabla\cdot(\rho{\bf v}) = 0$

(7.1.1.10)¶ ${\partial (\rho{\bf v})\over\partial t} + \nabla \cdot (\rho {\bf v}{\bf v}^T) + \nabla p = 0$

(7.1.1.11)¶ ${\partial E\over\partial t} + \nabla\cdot\left({\bf v}\left(E + p \right)\right) = 0$

The momentum equation can be further simplified by expanding the parentheses and using the continuity equation:

${\partial (\rho{\bf v})\over\partial t} + \nabla \cdot (\rho {\bf v}{\bf v}^T) + \nabla p = 0 \underbrace{\left( {\partial \rho\over \partial t} + \nabla\cdot(\rho{\bf v})\right)}_0 {\bf v} + \rho\left({\partial {\bf v}\over\partial t} + {\bf v}\cdot \nabla {\bf v}\right) + \nabla p = 0$

(7.1.1.12)¶ $\rho\left({\partial {\bf v}\over\partial t} + {\bf v}\cdot \nabla{\bf v}\right) + \nabla p = 0$

Where we used:

$\left[\nabla \cdot (\rho {\bf v}{\bf v}^T)\right]^i = \partial_j (\rho v^iv^j) = v^i \partial_j (\rho v^j) + \rho v^j \partial_j v^i = \left[{\bf v}\nabla \cdot (\rho {\bf v}) + \rho {\bf v}\cdot\nabla{\bf v}\right]^i$

Alternative Derivation¶

We can also take the non-relativistic limit in the stress energy tensor:

$T^{00} \to \rho c^2 T^{0i} = T^{i0} \to {1\over c} \rho c^2 v^i T^{ij} \to \rho v^iv^j + p \delta^{ij}$

and plug it into the equation (7.1.1.1). For $\mu=0$ we get the equation of continuity:

$\partial_\nu T^{0\nu} = 0 \partial_0 T^{00} + \partial_i T^{0i} = 0 \partial_t\left({1\over c} \rho c^2\right) + \partial_i\left({1\over c}\rho c^2 v^i \right) = 0 \partial_t \rho + \partial_i\left(\rho v^i \right) = 0$

and for $\mu=i$ we get the momentum equation:

$\partial_\nu T^{i\nu} = 0 \partial_0 T^{i0} + \partial_j T^{ij} = 0 \partial_t \left({1\over c^2}\rho c^2 v^i\right) + \partial_j \left( \rho v^iv^j + p \delta^{ij} \right) = 0 \partial_t \left(\rho v^i\right) + \partial_j \left( \rho v^iv^j + p \delta^{ij} \right) = 0$

However, in order to derive the equation for energy $E$ , one needs to take into account the full relativistic stress energy tensor, see the previous section for details.

Energy Equation¶

The energy equation can also be derived from thermodynamic and the other two Euler equations. We have the following two Euler equations:

$\partial_t\rho + \partial_i(\rho u^i) = 0 \rho\partial_t u^i + \rho u^j\partial_j u^i + \delta^{ij}\partial_j p = 0$

We’ll need the following formulas:

$\partial_t (u_i u^i) = (\partial_t u_i) u^i + u_i \partial_t u^i = (\partial_t u_i)\delta^{ij} u_j + u_i \partial_t u^i = = (\partial_t u_i\delta^{ij}) u_j + u_i \partial_t u^i = (\partial_t u^j) u_j + u_i \partial_t u^i = 2 u_i \partial_t u^i \partial_j (u_i u^i) = 2 u_i \partial_j u^i \partial_t\rho =- \partial_i(\rho u^i) \partial_t u^i =- u^j\partial_j u^i - {\delta^{ij}\over\rho}\partial_j p - u^j\partial_j p + \partial_t(\rho U) = = - {\d p \over\d t} +\partial_t p + \partial_t(\rho U) = = - {\d p \over\d t} +\partial_t (\rho U + p) = = - {\d p \over\d t} +{\d\over\d t} (\rho U + p) -u^j\partial_j (\rho U + p)= = - {\d p \over\d t} +{\d\rho\over\d t} \left(U + {p\over\rho}\right) +\rho{\d\over\d t} \left(U + {p\over\rho}\right) -u^j\partial_j (\rho U + p)= = - {\d p \over\d t} +{\d\rho\over\d t} \left(U + {p\over\rho}\right) +\rho{\d\over\d t} \left(U + {p\over\rho}\right) + (\rho U + p)\partial_j u^j -\partial_j (\rho U u^j + p u^j) = = \left[\rho {\d\over\d t}\left(U + {p\over\rho}\right) - {\d p\over\d t} \right] + \left(U + {p\over\rho}\right)\left[ {\d\rho\over\d t} + \rho \partial_j u^j \right] -\partial_j (\rho U u^j + p u^j) = = - \partial_j(\rho U u^j + p u^j) 0 = \d Q = T\d S = \d U + p\d V = \d (U + pV) - V\d p = \d\left(U+{p\over\rho}\right) - {1\over \rho}\d p = \d H - {1\over \rho}\d p$

where $V = {1\over\rho}$ is the specific volume and $H = U+{p\over\rho}$ is entalphy (heat content).

Then:

$\partial_t E = = \partial_t (\half \rho u_i u^i + \rho U) = = \half u_i u^i \partial_t \rho +\half\rho\partial_t(u_i u^i) + \partial_t(\rho U) = = -\half u_i u^i \partial_j(\rho u^j) +\rho u_i\partial_t u^i + \partial_t(\rho U) = = -\half u_i u^i \partial_j(\rho u^j) - \rho u_i u^j\partial_j u^i - u_i\delta^{ij}\partial_j p + \partial_t(\rho U) = = -\half u_i u^i \partial_j(\rho u^j) - \half\rho u^j\partial_j (u_i u^i) - u_i\delta^{ij}\partial_j p + \partial_t(\rho U) = = -\half\partial_j(\rho u_i u^i u^j) - u^j\partial_j p + \partial_t(\rho U) = = -\half\partial_j(\rho u_i u^i u^j) - \partial_j(\rho U u^j + p u^j) = = -\partial_j\left(u^j\left(\half\rho u_i u^i+\rho U + p \right)\right) = = -\partial_j\left(u^j\left(E + p \right)\right)$

so:

$\partial_t E + \partial_j\left(u^j\left(E + p \right)\right) = 0 {\partial E\over\partial t} + \nabla\cdot\left({\bf u}\left(E + p \right)\right) = 0$

7.1.2. Navier-Stokes Equations¶

We start with the following nonrelativistic components of the stress energy tensor:

$T^{00} \to \rho c^2 T^{0i} = T^{i0} \to {1\over c} \rho c^2 v^i T^{ij} \to \rho v^iv^j - \sigma^{ij}$

where $\sigma^{ij} = -p \delta^{ij} + \mathds{T}$ (more below) and plug it into the equation (7.1.1.1). For $\mu=0$ we get the equation of continuity as for perfect fluids:

$\partial_\nu T^{0\nu} = 0 \partial_0 T^{00} + \partial_i T^{0i} = 0 \partial_t\left({1\over c} \rho c^2\right) + \partial_i\left({1\over c}\rho c^2 v^i \right) = 0 \partial_t \rho + \partial_i\left(\rho v^i \right) = 0$

and for $\mu=i$ we get the momentum equation:

$\partial_\nu T^{i\nu} = f^i \partial_0 T^{i0} + \partial_j T^{ij} = f^i \partial_t \left({1\over c^2}\rho c^2 v^i\right) + \partial_j \left( \rho v^iv^j - \sigma^{ij} \right) = f^i \partial_t \left(\rho v^i\right) + \partial_j \left( \rho v^iv^j - \sigma^{ij} \right) = f^i$

By using the continuity equation in the momentum equation (as in perfect fluids), we get:

$\rho\left(\partial_t v^i + v^i \partial_j v^j\right) - \partial_j\sigma^{ij} = f^i$

This is sometimes called the Cauchy momentum equation:

$\rho\left({\partial {\bf v}\over\partial t} +{\bf v}\cdot\nabla{\bf v} \right) = \nabla \cdot \mathds{\sigma} + {\bf f}$

where the stress tensor $\sigma$ can be written as:

$\sigma=-p\mathds{1} + \mathds{T}$

and we get the Navier-Stokes equations:

$\rho\left({\partial {\bf v}\over\partial t} +{\bf v}\cdot\nabla{\bf v} \right) = -\nabla p + \nabla \cdot \mathds{T} + {\bf f}$

Those are the most general equations. If we assume some more things about the fluid, they can be further simplified.

For Newtonian fluids, we want $\mathds{T}$ to be isotropic, linear in strain rates and it’s divergence zero for fluid at rest. It follows that the only way to write the tensor under these conditions is:

$T_{ij} = 2\mu\epsilon_{ij} + \delta_{ij} \lambda \nabla\cdot{\bf v}$

where the strain rate is:

$\epsilon_{ij}={1\over 2}\left(\partial_j v_i+\partial_i v_j\right)$

The trace of $\mathds{T}$ is:

$\Tr \mathds{T} = T_{ii} = 2\mu\epsilon_{ii} + \delta_{ii} \lambda \nabla\cdot{\bf v} = (2\mu + 3 \lambda) \nabla\cdot{\bf v}$

Note that $\mathds{T}$ has zero trace, which is automatically satisfied for incompressible flow ( $\nabla\cdot{\bf v}=0$ ), but for compressible flow this imposes:

$\lambda = -{2\over 3}\mu$

The divergence of the tensor is:

$\partial_j T_{ij} =2\mu\partial_j\epsilon_{ij} + \partial_j\delta_{ij} \lambda \nabla\cdot{\bf v} =\mu\partial_j\partial_j v_i+\mu\partial_i \nabla\cdot{\bf v} + \lambda \partial_i \nabla\cdot{\bf v} =\mu\partial_j\partial_j v_i+(\mu+\lambda)\partial_i \nabla\cdot{\bf v}$

or in vector form (these are usually called the compressible Navier-Stokes equations):

$\nabla \cdot \mathds{T} =\mu\nabla^2{\bf v}+(\mu+\lambda)\nabla \nabla\cdot{\bf v}$

For incompressible fluid we have $\nabla\cdot\bf v=0$ , so we get the incompressible Navier-Stokes equations:

$\nabla \cdot \mathds{T} =\mu\nabla^2{\bf v}$

and for a perfect fluid we have no viscosity, e.g. $\mu=0$ , then we get the Euler equations (for perfect fluid):

$\nabla \cdot \mathds{T}=0$

7.1.3. Incompressible Equations¶

Incompressible flow means that the material derivative of density is zero:

(7.1.3.1)¶ ${\d \rho\over\d t} = {\partial \rho\over \partial t} + {\bf v}\cdot\nabla\rho = 0\,.$

Putting this into the equation of continuity (7.1.1.9) one obtains $\rho\nabla\cdot{\bf v}=0$ or equivalently:

(7.1.3.2)¶ $\nabla\cdot{\bf v}=0\,.$

But also (7.1.3.2) implies (7.1.3.1), so these two equations are equivalent: the divergence of the velocity field is zero if and only if the material derivative of the density is zero.

Using the condition $\nabla\cdot{\bf v}=0$ in (7.1.1.9) and (7.1.1.12) we obtain:

$\nabla\cdot{\bf v}=0\,, {\partial \rho\over \partial t} + {\bf v}\cdot\nabla\rho = 0\,, \rho\left({\partial {\bf v}\over\partial t} + {\bf v}\cdot \nabla{\bf v}\right) + \nabla p = \mu\nabla^2{\bf v}\,.$

In addition to incompressibility, we can also assume a constant density $\rho(x, y, z)=\rho_0$ , then we obtain the incompressible Navier-Stokes equations:

(7.1.3.3)¶ $\nabla\cdot{\bf v}=0\,,$

(7.1.3.4)¶ $\rho_0\left({\partial {\bf v}\over\partial t} + {\bf v}\cdot \nabla{\bf v}\right) + \nabla p = \mu\nabla^2{\bf v}\,.$

For $\mu=0$ they become the incompressible Euler equations. At the given time step with known ${\bf v}$ and $p$ , the equation (7.1.3.4) is solved for ${\bf v}$ at the new time step. Then we solve for new $p$ as follows. Apply divergence to (7.1.3.4):

$\rho_0\nabla\cdot\left({\partial {\bf v}\over\partial t} + {\bf v}\cdot \nabla{\bf v}\right) + \nabla\cdot\nabla p = \mu\nabla\cdot\nabla^2{\bf v}\,, \rho_0\left({\partial (\nabla\cdot{\bf v})\over\partial t} + \nabla\cdot({\bf v}\cdot \nabla{\bf v})\right) + \nabla^2 p = \mu\nabla\cdot\nabla^2{\bf v}\,,$

now we use the following identities:

$\nabla\cdot({\bf v}\cdot \nabla{\bf v}) = \partial_i(v^j \partial_j v^i) = (\partial_i v^j) (\partial_j v^i) + v^j \partial_j \partial_i v^i = \Tr (\nabla {\bf v})^2 + {\bf v}\cdot\nabla(\nabla\cdot{\bf v})\,, \nabla\cdot\nabla^2{\bf v} = \partial_i\partial^j\partial_j v^i = \partial^j\partial_j \partial_i v^i = \nabla^2(\nabla\cdot {\bf v})\,,$

to get:

$\rho_0\left({\partial (\nabla\cdot{\bf v})\over\partial t} + \Tr (\nabla {\bf v})^2 + {\bf v}\cdot\nabla(\nabla\cdot{\bf v}) \right) + \nabla^2 p = \mu \nabla^2(\nabla\cdot {\bf v}) \,.$

Finally we use the equation (7.1.3.3) to simplify:

(7.1.3.5)¶ $-\nabla^2 p = \rho_0\Tr (\nabla {\bf v})^2\,,$

which is a Poisson equation for $p$ . Note again that $\Tr (\nabla {\bf v})^2 = (\partial_i v^j) (\partial_j v^i)$ . The equation (7.1.3.5) is then used to solve for $p$ at the new time step.

Divergence Free Velocity¶

Typically by propagating (7.1.3.4), we obtain a velocity ${\bf v}^*$ that is not divergence free. To make it so, we want to find such a divergence free ${\bf v}$ that is closest to ${\bf v}^*$ in the $L^2$ norm $\|{\bf v} - {\bf v}^*\| = \sqrt{\int ({\bf v}-{\bf v}^*)^2 \, \d^3 x}$ , in other words we want to find the $L^2$ projection onto the divergence free subspace, so we have to minimize the following functional:

$R[{\bf v}, \lambda] = \half \| {\bf v}-{\bf v}^* \|^2 - \int \lambda \nabla\cdot{\bf v}\, \d^3 x = \int \half({\bf v}-{\bf v}^*)^2 - \lambda \nabla\cdot{\bf v}\, \d^3 x\,,$

where we used a Langrange multiplier $\lambda=\lambda({\bf x})$ in the second term to impose the zero divergence on ${\bf v}={\bf v}({\bf x})$ for all points ${\bf x}$ (that is why $\lambda$ is a function of ${\bf x}$ and not a constant) and in the first term we ensure that ${\bf v}$ is as close as possible to the original field ${\bf v}^*$ in the $L^2$ sense. Let’s calculate the variation:

$\delta R[{\bf v}, \lambda] = \int ({\bf v}-{\bf v}^*)\cdot\delta {\bf v} - \lambda \nabla\cdot \delta{\bf v} - (\nabla\cdot{\bf v}) \delta\lambda \, \d^3 x = = \int ({\bf v}-{\bf v}^*)\cdot\delta {\bf v} + (\nabla\lambda) \cdot \delta{\bf v} - (\nabla\cdot{\bf v}) \delta\lambda \, \d^3 x +\int \lambda \delta{\bf v}\cdot{\bf n}\, \d S = = \int ({\bf v}-{\bf v}^* + \nabla\lambda)\cdot\delta {\bf v} - (\nabla\cdot{\bf v}) \delta\lambda \, \d^3 x +\int \lambda \delta{\bf v}\cdot{\bf n}\, \d S\,.$

From the condition $\delta R[{\bf v}, \lambda]=0$ and assuming the surface integral vanishes (i.e. either $\lambda=0$ or $\delta{\bf v}\cdot{\bf n}=0$ everywhere on the boundary) we obtain the two Euler-Lagrange equations:

(7.1.3.6)¶ ${\delta R[{\bf v}, \lambda] \over \delta {\bf v}} = {\bf v}-{\bf v}^* + \nabla\lambda = 0\,,$

(7.1.3.7)¶ ${\delta R[{\bf v}, \lambda] \over \delta \lambda} = -\nabla\cdot{\bf v} = 0\,.$

Applying divergence to (7.1.3.6) and using (7.1.3.7) we obtain:

(7.1.3.8)¶ $\nabla^2\lambda = \nabla\cdot{\bf v}^*\,.$

After solving this Poisson equation for $\lambda$ we can calculate the divergence free ${\bf v}$ from (7.1.3.6):

(7.1.3.9)¶ ${\bf v} = {\bf v}^* - \nabla\lambda\,.$

Time Discretization¶

The incompressible Euler equations are:

$\nabla\cdot{\bf v}=0\,,$

(7.1.3.10)¶ ${\partial {\bf v}\over\partial t} = -{\bf v}\cdot \nabla{\bf v} -{1\over\rho_0} \nabla p\,.$

We use first order time discretization:

(7.1.3.11)¶ $\nabla\cdot{\bf v}^{n}=0\,,$

(7.1.3.12)¶ $\nabla\cdot{\bf v}^{n+1}=0\,,$

(7.1.3.13)¶ ${{\bf v}^{n+1} - {\bf v}^n\over\Delta t} = -{\bf v}^n\cdot \nabla{\bf v}^n -{1\over\rho_0} \nabla p^{n+1}\,.$

The velocity at time steps $n$ and $n+1$ must be divergence free, per (7.1.3.11) and (7.1.3.12). The simplest discretization of (7.1.3.10) is to use an explicit scheme, so we evaluate the term ${\bf v}\cdot \nabla{\bf v}$ at the time step $n$ . Regarding the pressure term $\nabla p$ , if we evaluated it at the time step $n$ , then from (7.1.3.13) we could calculate ${\bf v}^{n+1}$ that would not be divergence free, per (7.1.3.12). So we are led to evaluate the pressure term at the time step $n+1$ , then all the equations (7.1.3.11), (7.1.3.12) and (7.1.3.13) can be satisfied.

To solve this system of equations, we use an operator splitting on (7.1.3.13), the most natural is probably the following:

(7.1.3.14)¶ ${{\bf v}^* - {\bf v}^n\over\Delta t} = -{\bf v}^n\cdot \nabla{\bf v}^n -{1\over\rho_0} \nabla p^n\,,$

(7.1.3.15)¶ ${{\bf v}^{n+1} - {\bf v}^*\over\Delta t} = -{1\over\rho_0} \nabla (p^{n+1}-p^n)\,.$

The first equation (7.1.3.14) is just like (7.1.3.13), except that the pressure term is evaluated at the time step $n$ , which forces us to change ${\bf v}^{n+1}$ into ${\bf v}^*$ , which is not divergence free. The second equation (7.1.3.15) is then uniquely given by the condition that the sum of (7.1.3.14) and (7.1.3.15) is equal to (7.1.3.13).

The equation (7.1.3.15) is equivalent to (7.1.3.9), with $\lambda = {\Delta t \over\rho_0}(p^{n+1}-p^n)$ , so this is an $L^2$ projection of ${\bf v}^*$ onto the divergence free subspace to obtain ${\bf v}^{n+1}$ , also sometimes called a pressure projection. We use the same method as was used to obtain the Poisson equation (7.1.3.8) for $\lambda$ , i.e. take a divergence and rearrange:

(7.1.3.16)¶ $\nabla^2 \lambda = \nabla^2\left({\Delta t \over\rho_0}(p^{n+1}-p^n)\right) = \nabla\cdot{\bf v}^*\,.$

One solves (7.1.3.14) for ${\bf v}^*$ , then the Poisson equation (7.1.3.16) for $\lambda$ (i.e. the pressure update $p^{n+1}-p^n$ ), and then one computes ${\bf v}^{n+1}$ using (7.1.3.15) (or equivalently (7.1.3.9)).

These equations are derived from Euler equations (7.1.3.11) and (7.1.3.12) using a time discretization and an operator splitting technique. The theory of the $L^2$ projection onto the divergence free subspace is not needed to derive these equations, but it helps with understanding of what is going on.

Note 1: the operator splitting of (7.1.3.13) into (7.1.3.14) and (7.1.3.15) is not unique. Another option is:

(7.1.3.17)¶ ${{\bf v}^* - {\bf v}^n\over\Delta t} = -{\bf v}^n\cdot \nabla{\bf v}^n \,,$

(7.1.3.18)¶ ${{\bf v}^{n+1} - {\bf v}^*\over\Delta t} = -{1\over\rho_0} \nabla p^{n+1}\,.$

The sum of (7.1.3.17) and (7.1.3.18) is still (7.1.3.13) and the equation (7.1.3.18) is still equivalent to (7.1.3.9), only this time with $\lambda = {\Delta t \over\rho_0}p^{n+1}$ . The Poisson equation then becomes:

(7.1.3.19)¶ $\nabla^2 \lambda = \nabla^2\left({\Delta t \over\rho_0}p^{n+1}\right) = \nabla\cdot{\bf v}^*\,.$

The only difference to the previous scheme is that now the $L^2$ norm of $\| {\bf v}^{n+1} - {\bf v}^* \| = \| \nabla\lambda \|$ is larger, because $\lambda$ now depends on the full pressure instead of the pressure difference, so ${\bf v}^*$ is not as close to ${\bf v}^{n+1}$ as in the previous scheme.

Note 2: By applying divergence to (7.1.3.17) we obtain:

${\nabla\cdot{\bf v}^*\over\Delta t} = -\nabla\cdot({\bf v}^n\cdot \nabla{\bf v}^n) = -\Tr (\nabla {\bf v}^n)^2\,,$

and substituting into (7.1.3.19) we obtain:

$\nabla^2\left({\Delta t \over\rho_0}p^{n+1}\right) = \nabla\cdot{\bf v}^* = -\Delta t\Tr (\nabla {\bf v}^n)^2\,,$

or

(7.1.3.20)¶ $-\nabla^2 p^{n+1} = \rho_0 \Tr (\nabla {\bf v}^n)^2\,,$

which is the discrete analog of the equation (7.1.3.5). The same result is obtained by applying a divergence to (7.1.3.14) and substituting into (7.1.3.16):

$\nabla^2\left({\Delta t \over\rho_0}(p^{n+1}-p^n)\right) = \nabla\cdot{\bf v}^* = -\Delta t\Tr (\nabla {\bf v}^n)^2-{\Delta t\over\rho_0}\nabla^2 p^n$

Which simplifies to (7.1.3.20).

7.1.4. Bernoulli’s Principle¶

Bernoulli’s principle works for a perfect fluid, so we take the Euler equations:

$\rho\left({\partial {\bf v}\over\partial t} +{\bf v}\cdot\nabla{\bf v} \right) = -\nabla p + {\bf f}$

and put it into a vertical gravitational field ${\bf f} = (0, 0, -\rho g)=-\rho g\nabla z$ , so:

$\rho\left({\partial {\bf v}\over\partial t} +{\bf v}\cdot\nabla{\bf v} \right) = -\nabla p - \rho g\nabla z\,,$

we divide by $\rho$ :

${\partial {\bf v}\over\partial t} +{\bf v}\cdot\nabla{\bf v} = -\nabla \left({p\over\rho} + g z\right)$

and use the identity ${\bf v}\cdot\nabla{\bf v}={1\over 2}\nabla v^2 + (\nabla \times {\bf v})\times{\bf v}$ :

${\partial {\bf v}\over\partial t} +{1\over 2}\nabla v^2+(\nabla \times {\bf v})\times{\bf v} +\nabla \left({p\over\rho} + g z\right)=0\,,$

so:

${\partial {\bf v}\over\partial t} +(\nabla \times {\bf v})\times{\bf v} +\nabla \left({v^2\over 2} + gz + {p\over\rho} \right)=0\,.$

If the fluid is moving, we integrate this along a streamline from the point $A$ to $B$ :

$\int {\partial {\bf v}\over\partial t} \cdot \d {\bf l} +\left[{v^2\over 2} + gz + {p\over\rho} \right]_A^B=0\,.$

So far we didn’t do any approximation (besides having a perfect fluid in a vertical gravitation field). Now we assume a steady flow, so ${\partial {\bf v}\over\partial t}=0$ and since points $A$ and $B$ are arbitrary, we get:

${v^2\over 2} + gz + {p\over\rho}={\rm const.}$

along the streamline. This is called the Bernoulli’s principle. If the fluid is not moving, we set ${\bf v}=0$ in the equations above and immediately get:

$gz + {p\over\rho}={\rm const.}$

The last equation then holds everywhere in the (nonmoving) fluid (as opposed to the previous equation that only holds along the streamline).

Hydrostatic Pressure¶

Let $p_1$ be the pressure on the water surface and $p_2$ the pressure $h$ meters below the surface. From the Bernoulli’s principle:

${p_1\over\rho} = g\cdot (-h) + {p_2\over \rho}$

so

$p_1 + h\rho g = p_2$

and we can see, that the pressure $h$ meters below the surface is $h\rho g$ plus the (atmospheric) pressure $p_1$ on the surface.

Torricelli’s Law¶

We want to find the speed $v$ of the water flowing out of the tank (of the height $h$ ) through a small hole at the bottom. The (atmospheric) pressure at the water surface and also near the small hole is $p_1$ . From the Bernoulli’s principle:

${p_1\over\rho} = {v^2\over 2} + g\cdot (-h) + {p_1\over \rho}$

so:

$v=\sqrt{2g h}$

This is called the Torricelli’s law.

Venturi Effect¶

A pipe with a cross section $A_1$ , pressure $p_1$ and the speed of a perfect liquid $v_1$ changes it’s cross section to $A_2$ , so the pressure changes to $p_2$ and the speed to $v_2$ . Given $\Delta p = p_1-p_2$ , $A_1$ and $A_2$ , calculate $v_1$ and $v_2$ .

We use the continuity equation:

$A_1 v_1 = A_2 v_2$

and the Bernoulli’s principle:

${v_1^2\over 2} + {p_1\over\rho} = {v_2^2\over 2} + {p_2\over\rho}$

so we have two equations for two unknowns $v_1$ and $v_2$ , after solving it we get:

$v_1 = A_2\sqrt{2\Delta p\over \rho(A_1^2-A_2^2)}$

$v_2 = A_1\sqrt{2\Delta p\over \rho(A_1^2-A_2^2)}$

Hagen-Poiseuille Law¶

We assume incompressible (but viscuous) Newtonean fluid (in no external force field):

$\rho\left({\partial {\bf v}\over\partial t} +{\bf v}\cdot\nabla{\bf v} \right) = -\nabla p + \mu\nabla^2{\bf v}$

flowing in the vertical pipe of radius $R$ and we further assume steady flow ${\partial {\bf v}\over\partial t}=0$ , axis symmetry $v_r=v_\theta=\partial_\theta(\cdots)=0$ and a fully developed flow $\partial_z v_z=0$ . We write the Navier-Stokes equations above in the cylindrical coordinates and using the stated assumptions, the only nonzero equations are:

$0=-\partial_r p$

$0=-\partial_z p+\mu{1\over r}\partial_r(r\partial_r v_z)$

from the first one we can see the $p=p(z)$ is a function of $z$ only and we can solve the second one for $v_z=v_z(r)$ :

$v_z(r) = {1\over 4\mu}(\partial_z p) r^2 + C_1\log r + C_2$

We want $v_z(r=0)$ to be finite, so $C_1=0$ , next we assume the no slip boundary conditions $v_z(r=R)=0$ , so $C_2 = -{1\over 4\mu}(\partial_z p) R^2$ and we get the parabolic velocity profile: