Leonhard Euler

Figure 1: Leonhard Euler (1707-1783), source: Wiki)
  • Inventor and believer of absolute continuum hypotheses

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

  • Impact on today’s formalization convention

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

  • Groundbreaking work in mathematical modeling

Euler equations for fluid mechanics of frictionless elastic fluid flow. – Leonhard Euler: 1757

Motivation

Modeling the flow around an obstacle - here an airfoil:

Figure 2: Idealized flow overlay, Source: tec-science)

Motivation

Modeling the flow around an obstacle - here an airfoil:


Figure 3: Hele-Shaw flow past an inclined airfoil, source: An Album of Fluid Motion)

1 Incompressible Euler

Starting point

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!

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.

Strategy

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!

Incompressibility condition

Consider a medium (i) and physical situation (ii) subject to the following two assumptions:

  1. \(|\frac{1}{\rho} \frac{D}{Dt} \rho | << |\nabla \cdot \mathbf v |\)

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

Incompressibility condition

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

Incompressibility condition

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.

Ideal fluid approximation

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

Ideal fluid approximation

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.

2 Decomposition of the velocity gradient

3 Velocity gradient

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

Velocity gradient

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

Decomposition of the velocity field

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

4 Vorticity dynamics

Definition of the vorticity

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?

Lamb representation

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

Vorticity dynamics

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

Irrotational motion

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!

Potential flow

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

Outlook

Potential flow does not capture separation behavior!


Figure 4: Flow around airfoil with separation)

5 Bernoulli equation

Daniel Bernoulli

Figure 5: Daniel Bernoulli (1700 - 1782), source: Wiki
  • Child of an academic family

Father, brother and uncle were also mathematicians. Daniel, however, started off by studying medicine.

  • Gamma function

Bernoulli proposed an expression for the Gamma function

  • Riccati equation

Daniel Bernoulli became famous with his work on Riccatti equations.

  • Decomposition into translationary and rotational motion

  • many more contributions …

Bernoulli integral

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

Bernoulli integral

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

Bernoulli integral

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

Bernoulli integral

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

Bernoulli integral

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

Bernoulli integral

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.

Bernoulli equation

All in all: For

  • a medium in a stationary regime, hence \(\partial_t \mathbf v = 0\),
  • a barotropic medium for which \(P=\tfrac{p}{\rho}\), e.g. water, and
  • mass forces to be given as gravitational acceleration, hence \(U=gz\)

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

Bernoulli equation in a rotating reference frame

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

  • stationary (\(\partial_t \mathbf v =0\)),
  • irrotational flow (\(\nabla \times \mathbf{v} = 0\)),
  • and a fluid at rest w.r.t. the rotating reference frame:

\[ \frac{p}{\rho} - \frac{\left|\mathbf{\Omega}\right|^2\left|\mathbf{x}\right|^2}{2} + gz = \text{const}. \]

Bernoulli equation in a rotating reference frame

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.