.. index:: fluid dynamics Fluid Dynamics ============== .. index:: pair: stress-energy; tensor 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: .. math:: {\d p^\alpha\over\d V} = -T^\alpha{}_\beta u^\beta and the energy density that he measures is: .. math:: \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 .. math:: { \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}$: .. math:: {\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 .. math:: \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: .. math:: \partial_\nu T_\mu{}^\nu + f_\mu = 0 where .. math:: 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: .. math:: \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: .. math:: {\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 ~~~~~~~~~~~~~~ Perfect fluids have no heat conduction ($T^{i0} = T^{0i} = 0$) and no viscosity ($T^{ij} = p\one$), so in the comoving frame: .. math:: 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: .. math:: 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: .. math:: :label: conv1 \partial_\nu T^{\mu\nu} = 0 .. math:: :label: conv2 \partial_\mu(nu^\mu) = 0 The equation :eq:`conv2` gives: .. math:: \partial_t (n\gamma) + \partial_i(n v^i \gamma) = 0 .. math:: :label: continuity-relat \partial_t (n m\gamma) + \partial_i(n m v^i \gamma) = 0 .. math:: :label: continuity1 \partial_t (n m c^2\gamma) + \partial_i(n m c^2 v^i \gamma) = 0 The equation :eq:`conv1` gives for $\mu=0$: .. math:: \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 .. math:: :label: energy-relat2 \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 :eq:`continuity1` from :eq:`energy-relat2`: .. math:: \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: .. math:: 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 :eq:`energy-relat2`: .. math:: :label: energy-relat \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 :eq:`continuity-relat`: .. math:: :label: continuity-relat3 \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 :eq:`conv1` for $\mu=i$ gives: .. math:: \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 .. math:: :label: momentum-relat \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 :eq:`continuity-relat3`, :eq:`momentum-relat` and :eq:`energy-relat` are the correct relativistic equations for the perfect fluid (no approximations were done). We can take either :eq:`continuity-relat3` or :eq:`energy-relat2` 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$): .. math:: \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: .. math:: :label: euler_continuity {\partial \rho\over \partial t} + \nabla\cdot(\rho{\bf v}) = 0 .. math:: :label: euler_momentum {\partial (\rho{\bf v})\over\partial t} + \nabla \cdot (\rho {\bf v}{\bf v}^T) + \nabla p = 0 .. math:: :label: euler_energy {\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: .. math:: {\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 .. math:: :label: euler_momentum2 \rho\left({\partial {\bf v}\over\partial t} + {\bf v}\cdot \nabla{\bf v}\right) + \nabla p = 0 Where we used: .. math:: \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: .. math:: 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 :eq:`conv1`. For $\mu=0$ we get the equation of continuity: .. math:: \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: .. math:: \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: .. math:: \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: .. math:: \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: .. math:: \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: .. math:: \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 Navier-Stokes Equations ----------------------- We start with the following nonrelativistic components of the stress energy tensor: .. math:: 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 :eq:`conv1`. For $\mu=0$ we get the equation of continuity as for perfect fluids: .. math:: \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: .. math:: \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: .. math:: \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: .. math:: \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: .. math:: \sigma=-p\mathds{1} + \mathds{T} and we get the Navier-Stokes equations: .. math:: \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: .. math:: T_{ij} = 2\mu\epsilon_{ij} + \delta_{ij} \lambda \nabla\cdot{\bf v} where the strain rate is: .. math:: \epsilon_{ij}={1\over 2}\left(\partial_j v_i+\partial_i v_j\right) The trace of $\mathds{T}$ is: .. math:: \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: .. math:: \lambda = -{2\over 3}\mu The divergence of the tensor is: .. math:: \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): .. math:: \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: .. math:: \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): .. math:: \nabla \cdot \mathds{T}=0 Incompressible Equations ------------------------ Incompressible flow means that the material derivative of density is zero: .. math:: :label: mat_deriv_zero {\d \rho\over\d t} = {\partial \rho\over \partial t} + {\bf v}\cdot\nabla\rho = 0\,. Putting this into the equation of continuity :eq:`euler_continuity` one obtains $\rho\nabla\cdot{\bf v}=0$ or equivalently: .. math:: :label: div_free \nabla\cdot{\bf v}=0\,. But also :eq:`div_free` implies :eq:`mat_deriv_zero`, 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 :eq:`euler_continuity` and :eq:`euler_momentum2` we obtain: .. math:: \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: .. math:: :label: incomp_euler_1 \nabla\cdot{\bf v}=0\,, .. math:: :label: incomp_euler_2 \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 :eq:`incomp_euler_2` is solved for ${\bf v}$ at the new time step. Then we solve for new $p$ as follows. Apply divergence to :eq:`incomp_euler_2`: .. math:: \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: .. math:: \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: .. math:: \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 :eq:`incomp_euler_1` to simplify: .. math:: :label: incomp_euler_3 -\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 :eq:`incomp_euler_3` is then used to solve for $p$ at the new time step. Divergence Free Velocity ~~~~~~~~~~~~~~~~~~~~~~~~ Typically by propagating :eq:`incomp_euler_2`, 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: .. math:: 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: .. math:: \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: .. math:: :label: lambda_eq1a {\delta R[{\bf v}, \lambda] \over \delta {\bf v}} = {\bf v}-{\bf v}^* + \nabla\lambda = 0\,, .. math:: :label: lambda_eq1b {\delta R[{\bf v}, \lambda] \over \delta \lambda} = -\nabla\cdot{\bf v} = 0\,. Applying divergence to :eq:`lambda_eq1a` and using :eq:`lambda_eq1b` we obtain: .. math:: :label: lambda_eq2 \nabla^2\lambda = \nabla\cdot{\bf v}^*\,. After solving this Poisson equation for $\lambda$ we can calculate the divergence free ${\bf v}$ from :eq:`lambda_eq1a`: .. math:: :label: lambda_eq3 {\bf v} = {\bf v}^* - \nabla\lambda\,. Time Discretization ~~~~~~~~~~~~~~~~~~~ The incompressible Euler equations are: .. math:: \nabla\cdot{\bf v}=0\,, .. math:: :label: eul01 {\partial {\bf v}\over\partial t} = -{\bf v}\cdot \nabla{\bf v} -{1\over\rho_0} \nabla p\,. We use first order time discretization: .. math:: :label: eul11 \nabla\cdot{\bf v}^{n}=0\,, .. math:: :label: eul12 \nabla\cdot{\bf v}^{n+1}=0\,, .. math:: :label: eul13 {{\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 :eq:`eul11` and :eq:`eul12`. The simplest discretization of :eq:`eul01` 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 :eq:`eul13` we could calculate ${\bf v}^{n+1}$ that would not be divergence free, per :eq:`eul12`. So we are led to evaluate the pressure term at the time step $n+1$, then all the equations :eq:`eul11`, :eq:`eul12` and :eq:`eul13` can be satisfied. To solve this system of equations, we use an operator splitting on :eq:`eul13`, the most natural is probably the following: .. math:: :label: eul21 {{\bf v}^* - {\bf v}^n\over\Delta t} = -{\bf v}^n\cdot \nabla{\bf v}^n -{1\over\rho_0} \nabla p^n\,, .. math:: :label: eul22 {{\bf v}^{n+1} - {\bf v}^*\over\Delta t} = -{1\over\rho_0} \nabla (p^{n+1}-p^n)\,. The first equation :eq:`eul21` is just like :eq:`eul13`, 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 :eq:`eul22` is then uniquely given by the condition that the sum of :eq:`eul21` and :eq:`eul22` is equal to :eq:`eul13`. The equation :eq:`eul22` is equivalent to :eq:`lambda_eq3`, 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 :eq:`lambda_eq2` for $\lambda$, i.e. take a divergence and rearrange: .. math:: :label: eul23 \nabla^2 \lambda = \nabla^2\left({\Delta t \over\rho_0}(p^{n+1}-p^n)\right) = \nabla\cdot{\bf v}^*\,. One solves :eq:`eul21` for ${\bf v}^*$, then the Poisson equation :eq:`eul23` for $\lambda$ (i.e. the pressure update $p^{n+1}-p^n$), and then one computes ${\bf v}^{n+1}$ using :eq:`eul22` (or equivalently :eq:`lambda_eq3`). These equations are derived from Euler equations :eq:`eul11` and :eq:`eul12` 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 :eq:`eul13` into :eq:`eul21` and :eq:`eul22` is not unique. Another option is: .. math:: :label: eul31 {{\bf v}^* - {\bf v}^n\over\Delta t} = -{\bf v}^n\cdot \nabla{\bf v}^n \,, .. math:: :label: eul32 {{\bf v}^{n+1} - {\bf v}^*\over\Delta t} = -{1\over\rho_0} \nabla p^{n+1}\,. The sum of :eq:`eul31` and :eq:`eul32` is still :eq:`eul13` and the equation :eq:`eul32` is still equivalent to :eq:`lambda_eq3`, only this time with $\lambda = {\Delta t \over\rho_0}p^{n+1}$. The Poisson equation then becomes: .. math:: :label: eul33 \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 :eq:`eul31` we obtain: .. math:: {\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 :eq:`eul33` we obtain: .. math:: \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 .. math:: :label: eul5 -\nabla^2 p^{n+1} = \rho_0 \Tr (\nabla {\bf v}^n)^2\,, which is the discrete analog of the equation :eq:`incomp_euler_3`. The same result is obtained by applying a divergence to :eq:`eul21` and substituting into :eq:`eul23`: .. math:: \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 :eq:`eul5`. Bernoulli's Principle --------------------- Bernoulli's principle works for a perfect fluid, so we take the Euler equations: .. math:: \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: .. math:: \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$: .. math:: {\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}$: .. math:: {\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: .. math:: {\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$: .. math:: \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: .. math:: {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: .. math:: 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: .. math:: {p_1\over\rho} = g\cdot (-h) + {p_2\over \rho} so .. math:: 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: .. math:: {p_1\over\rho} = {v^2\over 2} + g\cdot (-h) + {p_1\over \rho} so: .. math:: 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: .. math:: A_1 v_1 = A_2 v_2 and the Bernoulli's principle: .. math:: {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: .. math:: v_1 = A_2\sqrt{2\Delta p\over \rho(A_1^2-A_2^2)} .. math:: v_2 = A_1\sqrt{2\Delta p\over \rho(A_1^2-A_2^2)} .. index:: pair: Hagen-Poiseuille; Law Hagen-Poiseuille Law ~~~~~~~~~~~~~~~~~~~~ We assume incompressible (but viscuous) Newtonean fluid (in no external force field): .. math:: \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: .. math:: 0=-\partial_r p .. math:: 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)$: .. math:: 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: .. math:: 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: .. math:: v_z(r) = {\Delta P\over 4\mu L}(R^2-r^2) We can now calculate the volumetric flow rate: .. math:: 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 = .. math:: ={\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.