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 \L\over\partial (\partial_\nu \eta_\rho)}

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 \L\over\partial (\partial_\nu \eta_\rho)}
    \partial_\mu \eta_\rho
    { \partial \L\over\partial (\partial_\nu \eta_\rho)}
    \partial_\nu\partial_\mu\eta_\rho + \partial_\mu\L

    { \partial \L\over\partial (\partial_\nu \eta_\rho)}
    \partial_\mu \eta_\rho
    + \partial_\mu\L


{ \partial \L\over\partial (\partial_\nu \eta_\rho)}
\partial_\mu \eta_\rho
-\L \delta_\mu{}^\nu
+ \partial_\mu\L

This can be written as:

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


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:

    = \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

    = \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:

(\partial_\nu T^{\mu\nu} = 0

(\partial_\mu(nu^\mu) = 0

The equation ( gives:

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

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

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

The equation ( 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

(\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 ( from (

\partial_t\left(\left(\rho c^2\gamma - n m c^2 + p {v^2\over c^2}
    \gamma\right) + \partial_i\left(\left(\rho c^2\gamma -n m c^2 + p
    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 (

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

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

(\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)
    \right) = 0

The equation ( 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

(\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 (, ( and ( are the correct relativistic equations for the perfect fluid (no approximations were done). We can take either ( or ( 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:

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

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

({\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

{\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

(\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 ( 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}
    \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).


\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)


\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.3. Incompressible Equations

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

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

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

(\nabla\cdot{\bf v}=0\,.

But also ( implies (, 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 ( and ( 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:

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

(\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 ( is solved for {\bf v} at the new time step. Then we solve for new p as follows. Apply divergence to (

\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})
    + \nabla^2 p = \mu \nabla^2(\nabla\cdot {\bf v}) \,.

Finally we use the equation ( to simplify:

(\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 ( is then used to solve for p at the new time step.

Divergence Free Velocity

Typically by propagating (, 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:

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

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

Applying divergence to ( and using ( we obtain:

(\nabla^2\lambda = \nabla\cdot{\bf v}^*\,.

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

({\bf v} = {\bf v}^* - \nabla\lambda\,.

Time Discretization

The incompressible Euler equations are:

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

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

We use first order time discretization:

(\nabla\cdot{\bf v}^{n}=0\,,

(\nabla\cdot{\bf v}^{n+1}=0\,,

({{\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 ( and ( The simplest discretization of ( 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 ( we could calculate {\bf v}^{n+1} that would not be divergence free, per ( So we are led to evaluate the pressure term at the time step n+1, then all the equations (, ( and ( can be satisfied.

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

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

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

The first equation ( is just like (, 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 ( is then uniquely given by the condition that the sum of ( and ( is equal to (

The equation ( is equivalent to (, 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 ( for \lambda, i.e. take a divergence and rearrange:

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

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

These equations are derived from Euler equations ( and ( 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 ( into ( and ( is not unique. Another option is:

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

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

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

(\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 ( 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 ( 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\,,


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

which is the discrete analog of the equation ( The same result is obtained by applying a divergence to ( and substituting into (

\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.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\,,


{\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}


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}


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:

v_z(r) = {1\over 4\mu}(-\partial_z p) (R^2-r^2)

Assuming that the pressure decreases linearly across the length of the pipe, we have -\partial_z p = {\Delta P\over L} and we get:

v_z(r) = {\Delta P\over 4\mu L}(R^2-r^2)

We can now calculate the volumetric flow rate:

Q = {\d V\over\d t} ={\d\over \d t}\int z\, \d S =\int {\d z\over \d t} \d S =\int v_z \,\d S =\int_0^{2\pi}\int_0^R v_z\, r\, \d r\,\d\phi =

={\Delta P\pi\over 2\mu L}\int_0^R (R^2-r^2) r\, \d r ={\Delta P \pi R^4\over 8 \mu L}

so we can see that it depends on the 4th power of R. This is called the Hagen-Poiseuille law.