Incompressible Euler equations
Methods for Model-Based Development in CE
2025-04-25
Ich lege also fest, dass inmitten der Bewegung des Fluids vom Fluid her kein freier Raum zurückgelassen wird, sondern die Kontinuität in ihm beständig erhalten bleibt. – Leonhard Euler: 1761
Inventor of \(f(x)\), \(e\), \(i\), for example: \(f(x) = (x+a)^n= \sum P(\{ (e_1,...,e_N) \}) = {N \choose k} \cdot p^kq^{N-k}\)
Euler equations for fluid mechanics of frictionless elastic fluid flow. – Leonhard Euler: 1757
Modeling the flow around an obstacle - here an airfoil:
Modeling the flow around an obstacle - here an airfoil:
In the previous section, we derived mass and momentum balance according to
\[ \begin{aligned} \partial_t \rho + \nabla \cdot \left( \rho \mathbf v \right) & = 0 \\ \partial_t \left( \rho \mathbf v \right) + \nabla \cdot \mathbf \Pi & = \rho \mathbf b, \end{aligned} \]
with * \(\rho \mathbf v\) denoting the flux density vector with components, and * \(\rho v_i\) and \(\mathbf \Pi\) denoting the second order tensor \(\rho \mathbf v \otimes \mathbf v - \mathbf \sigma\) with components \(\Pi_{ij} = \rho v_i v_j - \sigma_{ij}\).
Observation: We now have four equations for 13 unknown \(\rho, v_i, \sigma_{ij}\).
Note, that considering the mass balance alone, we have one equation for four unknowns: \(\rho, v_i\), which is also undetermined. Considering another balance law, here the momentum balance, however, doesn’t help to close the problem. We need a different strategy!
Let’s characterize the conditions that our models should be valid in by narrowing down material properties of the continuous medium , and / or the physical regime that we would like our mathematical process model to cover.
Such specialization strategies aim at enriching the mathematical model with additional closure relations. We hence loose generality of the model, yet gain a necessary prerequisite for well posedness, hence solvability.
Assume the Cauchy stress simplifies into \[ \mathbf \sigma = - p \mathbf I + \mathbf \tau, \]
in which
\(p\) stands for the pressure, hence stresses due to isotropic volume change, and
\(\tau := \sigma - i \mathbf I\) denotes the so-called deviatoric part of \(\mathbf \sigma\).
For now, we simply accept the fact that this decomposition makes sense. A detailed discussion of the Cauchy stress tensor and its decomposition will follow later.
With that decomposition, we can re-write mass and momentum balance as
\[ \begin{aligned} \tfrac{1}{\rho} \tfrac{D}{Dt} \rho + \nabla \cdot \mathbf v & = 0 \\ \tfrac{D}{Dt} \mathbf v & = - \tfrac{1}{\rho} \nabla p + \tfrac{1}{\rho} \nabla \cdot \mathbf \tau + \mathbf b. \end{aligned} \]
Conclusion: This does not yet help to close the system!
Consider a medium (i) and physical situation (ii) subject to the following two assumptions:
\(|\frac{1}{\rho} \frac{D}{Dt} \rho | << |\nabla \cdot \mathbf v |\)
\(|\nabla \cdot \mathbf \tau | << | \nabla p|\).
Assumption (i) means that the total change in density of the medium is much smaller than the dominant term in \(\nabla \cdot \mathbf v = \partial_x u + \partial_y v + \partial_z w\).
Assumption (ii) means that the deviatoric stresses can be neglected in the momentum flux density vector \(\mathbf \Pi\).
We can use this to simplify mass and momentum balance equations:
Mass balance turns into \[ \begin{aligned} 0 = \underbrace{\frac{1}{\rho} \frac{D}{Dt} \rho}_{\rightarrow 0} + \nabla \cdot \mathbf v = \nabla \cdot \mathbf v \end{aligned} \]
Observation: The velocity field is divergence free (no sinks and sources in the domain).
We can use this to simplify mass and momentum balance equations:
Momentum balance turns into
\[ \begin{aligned} \rho \partial_t \mathbf v + \underbrace{\mathbf v \partial_t \rho + \nabla \cdot ( \rho \mathbf v) \mathbf v}_{\mathbf v \left(\partial_t \rho + \nabla \cdot ( \rho \mathbf v) \right) = 0} + \nabla \mathbf v (\rho \mathbf v) &= - \nabla p + \underbrace{\nabla \cdot \mathbf \tau}_{\rightarrow 0} + \rho \mathbf b \\ \Leftrightarrow \qquad \partial_t \mathbf v + \mathbf v \cdot \nabla \mathbf v &= - \frac{1}{\rho}\nabla p + \mathbf b \end{aligned} \]
Observation: As the deviatoric stresses drop out of the momentum balance, we are losing dependency on unknowns.
The above assumptions (i) and (ii) are called the ideal, incompressible fluid approximation.
The corresponding mathematical model is refers to the ideal, incompressible fluid model.
An alternative name ist the incompressible Euler equations.
It reads concisely written:
\[ \begin{aligned} \nabla \cdot \mathbf v & = 0\\ \partial_t \mathbf v + \mathbf v \cdot \nabla \mathbf v &= - \frac{1}{\rho}\nabla p + \mathbf b \end{aligned} \]
Alternatively, we can write it in index notation:
\[ \begin{aligned} \partial_i v_i & = 0\\ \partial_t v_i + v_j \partial_j v_i &= - \tfrac{1}{\rho}\partial_i p + b_i \end{aligned} \]
The above assumptions (i) and (ii) are called the ideal, incompressible fluid approximation.
The corresponding mathematical model is refers to the ideal, incompressible fluid model.
An alternative name ist the incompressible Euler equations.
It reads concisely written:
\[ \begin{aligned} \nabla \cdot \mathbf v & = 0\\ \partial_t \mathbf v + \mathbf v \cdot \nabla \mathbf v &= - \frac{1}{\rho}\nabla p + \mathbf b \end{aligned} \]
Assuming we know the constant density \(\rho\) and the body force, for instance when given as the gravitational acceleration \(\mathbf b = \mathbf g\), the incompressible Euler equations comprise a set of four equations in the four unknowns pressure \(p\) and velocity \(\mathbf v\).
We are able to solve it, given physically meaningful initial and boundary conditions.
Let’s consider the velocity field around an arbitrary point \(\mathbf x\). Expanding the velocity field in a Taylor series yields
\[ \mathbf v \left( \mathbf x + \mathbf r \right) = \mathbf v \left( \mathbf x \right) + \left( \nabla \mathbf v \right) \mathbf r + \text{higher order terms} \]
Here, \(\nabla \mathbf v\) is a second order tensor with components \(\partial_j v_i\).
The velocity gradient can be decomposed into a symmetric and an anti-symmetric part:
\[ \nabla \mathbf v = \mathbf W + \mathbf D, \]
with
\[ \mathbf W = \tfrac{1}{2} \left( \nabla \mathbf v - \nabla \mathbf v^T \right) \quad \text{and} \quad \mathbf D = \tfrac{1}{2} \left( \nabla \mathbf v + \nabla \mathbf v^T \right), \]
In components:
\[ W_{ij} = \tfrac{1}{2} \left( \partial_j v_i - \partial_i v_j \right) \quad \text{and} \quad D_{ij} = \tfrac{1}{2} \left( \partial_j v_i + \partial_i v_j \right). \]
Physical interpretation:
\(\mathbf D\) is called the strain-rate tensor that measures the strain or distorsion of a medium. Note, that the strain-rate tensor is symmetric and can therefore be rotated into a coordinate system, in which it is represented as a diagonal matrix. The diagonal entries \(\partial_i v_i\) account for dilation or compression in the three principal directions.
\(\mathbf W\) on the other hand denotes the spin tensor and accounts for local rotational motion or spin. \(\mathbf W\) is skew-symmetric and has zero entries on the diagonal, hence a vanishing trace. We saw earlier that
\[ \mathbf W \cdot \mathbf r = \mathbf w \times \mathbf r \]
with \(\mathbf w\) being the axial vector of \(\mathbf W\) defined as
\[ \mathbf w = \tfrac{1}{2} \text{curl} \, \mathbf v = \tfrac{1}{2} \nabla \times \mathbf v. \]
Recall, that \(\mathbf w \times \mathbf r\) indicates a rigid body rotation of a point \(\mathbf r\) around rotation axis \(\mathbf w = w \mathbf n\).
All in all, we get
\[ \underbrace{\mathbf v \left( \mathbf x + \mathbf r \right)}_{\text{velocity field}} = \underbrace{\mathbf v \left( \mathbf x \right)}_{\text{translational motion}} + \underbrace{\mathbf w \times \mathbf r }_{\text{rotational motion}} + \underbrace{\mathbf D \cdot \mathbf r}_{\text{internal deformation}}. \]
which means that the velocity field in a region around a certain location is given by a superposition of
local distorsion and
local rotation
We now define the so-called vorticity
\[\mathbf \omega = \nabla \times \mathbf v,\]
which is the axial vector for \(2 \mathbf W\), or \(2 \mathbf w\).
Question: How does the vorticity evolve under the ideal, incompressible fluid assumption?
By making use of the vector operator identity
\[ \left( \mathbf v \cdot \nabla \right) \mathbf v = \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right) \]
we can write Euler’s momentum equation of the Euler equations in the so-called Lamb representation (1895):
\[ \partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right) = \mathbf b - \tfrac{1}{\rho} \nabla p. \]
Applying the curl operation \(\nabla \times\) to both sides of the equation yields
\[ \nabla \times \left( \partial_t \mathbf v \right) + \nabla \times \left(\nabla \left( \tfrac{\mathbf v^2}{2} \right)\right) - \nabla \times \left(\mathbf v \times \left( \nabla \times \mathbf v \right)\right) = \nabla \times \mathbf b - \tfrac{1}{\rho} \nabla \times \nabla p, \]
which simplifies as we identify \(\mathbf \omega = \nabla \times \mathbf v\) and consider that \(\nabla \times \nabla *\) vanishes:
\[ \partial_t \mathbf \omega - \nabla \times \left( \mathbf v \times \mathbf \omega \right) = \nabla \times \mathbf b. \]
Next, we exploit the vector operator identity
\[ \nabla \times \left( \mathbf a \times \mathbf b \right) = \left( \mathbf b \cdot \nabla \right) \mathbf a - \left( \mathbf a \cdot \nabla \right) \mathbf b - \mathbf a \left( \nabla \cdot \mathbf b \right) - \mathbf b \left( \nabla \cdot \mathbf a \right), \]
by identifying \(\mathbf a = \mathbf v\) and \(\mathbf b = \mathbf \omega\).
Considering \(\nabla \cdot \left(\nabla \times * \right) = 0\) as well as continuity (\(\nabla \cdot \mathbf v=0\)) finally yields
\[ \partial_t \mathbf \omega = \left( \mathbf \omega \cdot \nabla \right) \mathbf v - \left( \mathbf v \cdot \nabla \right) \mathbf \omega + \nabla \times \mathbf b. \]
Alternatively, written in terms of the material derivative:
\[ \tfrac{D}{Dt} \mathbf \omega = \left( \mathbf \omega \cdot \nabla \right) \mathbf v + \nabla \times \mathbf b. \]
This means that total change of the vorticity evolution is due to vorticity production by means of existing vorticity in a non-uniform velocity field, and due to the vortex components of the mass force field \(\mathbf b\).
If the mass force field constitutes a potential, hence \(\mathbf b = \nabla B\), we get
\[ \tfrac{D}{Dt} \mathbf \omega = \left( \mathbf \omega \cdot \nabla \right) \mathbf v. \]
An important consequence of this is that if the mass force is a potential and the initial vorticity is zero everywhere, we get
\[ \tfrac{D}{Dt} \mathbf \omega = 0 \]
always!
Under the previous condition (mass force is a potential, no initial vorticity), the velocity field is characterized by a zero vorticity, hence \(\nabla \times \mathbf v = 0\).
Due to mass conservation, the velocity field also has a vanishing divergence \(\nabla \cdot \mathbf v=0\).
In such a situation, the velocity itself is given as a potential, meaning there is a scalar function \(\phi\), referred to as the velocity potential, such that
\[ \mathbf v = \nabla \phi \]
In combination with the continuity equation, we get
\[ 0 = \nabla \cdot \mathbf v = \nabla \cdot \nabla \phi = \triangle \phi \]
Remark
The velocity potential \(\phi\) follows a (linear) Laplace equation. This can be solved in a straight forward way via standard finite element techniques. The pressure has to be recovered in a second step by solving the (non-linear) Euler equation for pressure based on the known velocity field.
Check out this great resource to get some hands-on experience with implementing potential flow: Aero Python
Potential flow does not capture separation behavior!
Father, brother and uncle were also mathematicians. Daniel, however, started off by studying medicine.
Bernoulli proposed an expression for the Gamma function
Daniel Bernoulli became famous with his work on Riccatti equations.
Decomposition into translationary and rotational motion
many more contributions …
We saw that the Euler equations can be put into Lamb representation.
\[ \partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right) = \mathbf b - \tfrac{1}{\rho} \nabla p. \]
Let’s integrate this expression along a streamline \(\mathbf x(t)\) starting from coordinate \(\mathbf x_1\) up to \(\mathbf x_2\), which for simplicity are denoted by 1 and 2:
\[ \int_1^2 \left(\partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right)\right) \cdot d \mathbf x = \int_1^2 \mathbf b \cdot d \mathbf x - \int_1^2 \tfrac{1}{\rho} \nabla p \cdot d \mathbf x. \]
If the mass force is once again given as a potential \(\mathbf b = - \nabla U\), we get
\[ \int_1^2 \nabla \left( \tfrac{\mathbf v^2}{2} \right)\cdot d \mathbf x - \int_1^2 \mathbf b \cdot d \mathbf x = \int_1^2 \nabla \left( \tfrac{\mathbf v^2}{2} + U \right)\cdot d \mathbf x \]
We saw that the Euler equations can be put into Lamb representation.
\[ \partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right) = \mathbf b - \tfrac{1}{\rho} \nabla p. \]
Let’s integrate this expression along a streamline \(\mathbf x(t)\) starting from coordinate \(\mathbf x_1\) up to \(\mathbf x_2\), which for simplicity are denoted by 1 and 2:
\[ \int_1^2 \left(\partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right)\right) \cdot d \mathbf x = \int_1^2 \mathbf b \cdot d \mathbf x - \int_1^2 \tfrac{1}{\rho} \nabla p \cdot d \mathbf x. \]
The integral of \(\nabla \left( \tfrac{\mathbf v^2}{2} + U \right)\) can be evaluated directly:
\[ \int_1^2 \nabla \left( \tfrac{\mathbf v^2}{2} + U \right)\cdot d \mathbf x = \left[ \tfrac{\mathbf v^2}{2} + U \right]_1^2 \]
We saw that the Euler equations can be put into Lamb representation.
\[ \partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right) = \mathbf b - \tfrac{1}{\rho} \nabla p. \]
Let’s integrate this expression along a streamline \(\mathbf x(t)\) starting from coordinate \(\mathbf x_1\) up to \(\mathbf x_2\), which for simplicity are denoted by 1 and 2:
\[ \int_1^2 \left(\partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right)\right) \cdot d \mathbf x = \int_1^2 \mathbf b \cdot d \mathbf x - \int_1^2 \tfrac{1}{\rho} \nabla p \cdot d \mathbf x. \]
By introducing the pressure increment along the integration path \(d p = \nabla p \cdot d \mathbf x\), we can transform the integral of the pressure gradient and get
\[ \int_1^2 \tfrac{1}{\rho} \nabla p \cdot d \mathbf x = \int_{p_1}^{p_2} \tfrac{1}{\rho} dp \]
We saw that the Euler equations can be put into Lamb representation.
\[ \partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right) = \mathbf b - \tfrac{1}{\rho} \nabla p. \]
Let’s integrate this expression along a streamline \(\mathbf x(t)\) starting from coordinate \(\mathbf x_1\) up to \(\mathbf x_2\), which for simplicity are denoted by 1 and 2:
\[ \int_1^2 \left(\partial_t \mathbf v + \nabla \left( \tfrac{\mathbf v^2}{2} \right) - \mathbf v \times \left( \nabla \times \mathbf v \right)\right) \cdot d \mathbf x = \int_1^2 \mathbf b \cdot d \mathbf x - \int_1^2 \tfrac{1}{\rho} \nabla p \cdot d \mathbf x. \]
Let’s furthermore assume a vanishing time dependency:
\[ \int_1^2 \left(\partial_t \mathbf v \right) \cdot d \mathbf x = 0 \]
We are then left with the following general formulation of the Bernoulli equation:
\[ \left[ \tfrac{\mathbf v^2}{2} + U \right]_1^2 + \int_{p_1}^{p_2} \tfrac{dp}{\rho} = \int_1^2 \left( \mathbf v \times \left( \nabla \times \mathbf v \right)\right) \cdot d \mathbf x \]
In a so-called barotropic fluid the pressure is a function of density alone (\(p=p(\rho)\)) and vice versa (\(\rho=\rho(p)\)). As a consequence, we can introduce a pressure function
\[P(p):= \int_{p_0}^p \tfrac{d \tilde p}{\rho(\tilde p)},\]
such that
\[ \int_{p_1}^{p_2} \tfrac{dp}{\rho} = P(p_2) - P(p_1) =: P_2 - P_1 \]
We are then left with the following general formulation of the Bernoulli equation:
\[ \left[ P - \tfrac{\mathbf v^2}{2} + U \right]_1^2 = \underbrace{\int_1^2 \left( \mathbf v \times \left( \nabla \times \mathbf v \right)\right) \cdot d \mathbf x}_{=\text{vanishes along streamlines}} = 0 \]
The right hand side, finally vanishes along streamlines.
Note that the right hand side also vanishes for irrotational flow, such the Bernoulli equation also holds for arbitrary locations in space.
All in all: For
we result at the famous Bernoulli equation, which reads
\[ \tfrac{p}{\rho} + \tfrac{v^2}{2} + g z = const \]
Often the static pressure \(p\) is combined with the pressure due to gravitational acceleration into the piezometric pressure \(p^* = p + \rho g z\), hence
\[ p^* + \rho\tfrac{v^2}{2} = const \]
Let’s consider a reference frame that rotates at \(\mathbf \Omega\). The Lamb representation is then
\[ \partial_t \mathbf v + \nabla \left( \frac{\mathbf{v}^2}{2} - \underbrace{\frac{\left|\mathbf{\Omega}\right|^2\left|\mathbf{x}\right|^2}{2}}_{=\text{centrifugal force}} \right) + \underbrace{2 \mathbf{\Omega} \times \mathbf{v}}_{=\text{Coriolis force}}- \mathbf{v} \times \nabla \times \mathbf{v} = - \frac{1}{\rho}\nabla p + \mathbf b. \]
Following similar arguments than before, we derive the Bernoulli equation for
\[ \frac{p}{\rho} - \frac{\left|\mathbf{\Omega}\right|^2\left|\mathbf{x}\right|^2}{2} + gz = \text{const}. \]
Balancing a location on the rotational axis (\(x_1=0, z_1=0\)) with a point away from the axis (x_2 , z_2 $) yields
\[ \frac{p_2}{\rho} - \frac{\left|\mathbf{\Omega}\right|^2 x_2^2}{2} + gz_2 = \frac{p_1}{\rho}, \]
with \(p_1=p_2=p_\mathrm{ambient}\), we get
\[ z_2 = \frac{\left|\mathbf{\Omega}\right|^2 x_2^2}{2 g}, \]
and so the surface is curved to form a rotational paraboloid.
Continuum Mechanical Modeling for Simulation Science - Incompressible Euler equations