Where Total Pressure Dropped 1,080 Pa, Dissipation Was Zero — When the Bernoulli Constant Crosses Streamlines

In 1738 a father backdated his son's book by six years#

Daniel Bernoulli published Hydrodynamica in 1738. On the title page he described himself as "son of Johann," a gesture meant to patch things up with his father. Johann answered by publishing nearly the same material separately as Hydraulica, leaning on his printer to stamp the year 1732 on it. Six years ahead of his son, on paper.

A century and a half later Horace Lamb supplied the theorem. Integrate the Euler equations and Bernoulli's result falls out. The two books the father and son fought over were two faces of one expression. But that "integrate" carries a condition attached, and that condition is exactly what makes total-pressure contours get misread in CFD post-processing.

This post measures the condition with a single Rankine vortex. The answer first: inside the rotating core the total pressure drops 1,080 Pa, and viscous dissipation is exactly zero. Outside, where the total pressure is perfectly flat, dissipation is not zero. The two come out backwards.

Move the probe in and out of the core yourself in the simulation below.

Drag probe r and watch where the green $p_0$ curve on the right goes flat. At the same time watch the square on the left of the canvas, and see when it holds its shape and when it skews. The place where those two observations disagree is the subject of this post.

Integrating the Euler equations gives Bernoulli — in which direction?#

For steady inviscid flow at constant density with no body force, the momentum equation reads like this.

$\mathbf{u}$ is velocity, $p$ is pressure, $\rho$ is density. Apply a vector identity to the left side.

$\boldsymbol{\omega} = \nabla\times\mathbf{u}$ is the vorticity. Combine the two and collect the total pressure $p_0 = p + \tfrac{1}{2}\rho|\mathbf{u}|^2$, and one line is left.

This is the Crocco form of the Euler equation (incompressible, isentropic). When the right side is nonzero, total pressure varies in space. Which direction it varies in is the point.

$\mathbf{u}\times\boldsymbol{\omega}$ is perpendicular to $\mathbf{u}$. So dotting both sides with $\mathbf{u}$ makes the right side vanish.

Follow a streamline and $p_0$ is always constant. Vorticity or no vorticity, it makes no difference. That is the exact scope of Bernoulli's theorem. For $p_0$ to be the same everywhere, on the other hand, you need $\mathbf{u}\times\boldsymbol{\omega} = 0$, which in practice means irrotational flow. The difference between the two feuding books was never here. The condition sat in the footnote Lamb added.

The Rankine vortex: a core that spins, an exterior that does not#

The Rankine vortex is the example where this condition switches on and off inside one flow. Inside radius $a$ the fluid turns like a solid body, and outside it is a free vortex.

$\Omega$ is the angular velocity of the core, $u_\theta$ the circumferential velocity. Pressure comes from integrating the radial momentum balance $\mathrm{d}p/\mathrm{d}r = \rho u_\theta^2 / r$. Outside, the amount the pressure falls and the amount the dynamic pressure rises cancel exactly, pinning $p_0$ at $p_\infty$. Inside the core nothing cancels.

At $r = 0$ the deficit is $\rho\Omega^2 a^2$, that is, the square of the core-edge speed times $\rho$. With $\rho = 1.2$, $\Omega = 60\,\mathrm{s^{-1}}$ and $a = 0.5\,\mathrm{m}$, $u_{\max} = 30$ m/s and the deficit is 1,080 Pa.

Measuring both sides at the same radius in Python#

Take $\nabla p_0$ by central difference and match it against $\rho\,u_\theta\omega_z$. The dissipation function $\phi = \mu(2S_{r\theta})^2$ goes on the same ruler.

import math
 
RHO, MU = 1.2, 1.8e-5      # kg/m^3, Pa*s
OMEGA, A = 60.0, 0.5       # 1/s, m   -> u_max = 30 m/s
P_INF = 101325.0           # Pa
 
def rankine_velocity(r):
    return OMEGA * r if r < A else OMEGA * A * A / r
 
def rankine_spin(r):
    return 2.0 * OMEGA if r < A else 0.0
 
def rankine_pressure(r):
    if r >= A:
        return P_INF - 0.5 * RHO * (OMEGA * A * A / r) ** 2
    return P_INF - RHO * OMEGA**2 * A**2 + 0.5 * RHO * OMEGA**2 * r**2
 
def total_head(r):
    return rankine_pressure(r) + 0.5 * RHO * rankine_velocity(r) ** 2
 
def shear_dissipation(r):
    s = 0.0 if r < A else -OMEGA * A * A / r**2   # 2*S_rtheta
    return MU * s * s
 
def crocco_residual(r, h=1e-6):
    lhs = (total_head(r + h) - total_head(r - h)) / (2 * h)   # d(p0)/dr
    rhs = RHO * rankine_velocity(r) * rankine_spin(r)         # rho*(u x omega)_r
    return lhs, rhs
 
print("  r[m]   u[m/s]   p[Pa]      p0[Pa]     dp0/dr    rho*u*w    phi[W/m^3]")
for r in (0.10, 0.25, 0.40, 0.60, 1.00):
    lhs, rhs = crocco_residual(r)
    print("%6.2f %8.2f %10.1f %10.1f %9.1f %10.1f %11.4f"
          % (r, rankine_velocity(r), rankine_pressure(r), total_head(r),
             lhs, rhs, shear_dissipation(r)))
 
def bernoulli_gap(r1, r2):
    return total_head(r2) - total_head(r1)
 
print()
print("core  p0(0.00) - p0(%.2f) = %8.1f Pa" % (A, total_head(0.0) - total_head(A)))
print("outer p0(%.2f) - p0(1.00) = %8.1f Pa" % (A, total_head(A) - total_head(1.0)))
print("cross-streamline gap r=0.10 -> 0.40 : %8.1f Pa" % bernoulli_gap(0.10, 0.40))
print("along-streamline  gap r=0.40 -> 0.40 : %8.1f Pa" % bernoulli_gap(0.40, 0.40))
print("dissipation at r=0.25 (core)  : %.4f W/m^3" % shear_dissipation(0.25))
print("dissipation at r=0.60 (outer) : %.4f W/m^3" % shear_dissipation(0.60))

  r[m]   u[m/s]   p[Pa]      p0[Pa]     dp0/dr    rho*u*w    phi[W/m^3]
  0.10     6.00   100266.6   100288.2     864.0      864.0      0.0000
  0.25    15.00   100380.0   100515.0    2160.0     2160.0      0.0000
  0.40    24.00   100590.6   100936.2    3456.0     3456.0      0.0000
  0.60    25.00   100950.0   101325.0       0.0        0.0      0.0313
  1.00    15.00   101190.0   101325.0       0.0        0.0      0.0040
 
core  p0(0.00) - p0(0.50) =  -1080.0 Pa
outer p0(0.50) - p0(1.00) =      0.0 Pa
cross-streamline gap r=0.10 -> 0.40 :    648.0 Pa
along-streamline  gap r=0.40 -> 0.40 :      0.0 Pa
dissipation at r=0.25 (core)  : 0.0000 W/m^3
dissipation at r=0.60 (outer) : 0.0313 W/m^3

Press match outlet spread and the outlet total-pressure span of the two channels matches, so only a particle whose color changes along its path still counts as loss.

Particles in the upper channel keep their color. Only the lower channel changes color as it passes the body. The 1,080 Pa core deficit in the output just above sits with the channel whose color does not change. Along a streamline that gap is zero.

Columns four and five agree at every radius. The Crocco relation holds down to the decimal. And inside the core the total-pressure difference between $r = 0.10$ and $r = 0.40$ is 648 Pa. Between two angles at the same radius it is zero. Bernoulli applied along a streamline is right; applied across one it misses 648 Pa.

The three slopes line up with radius. At 0.10 m the gradient is 864 Pa/m, at 0.25 m it is 2,160, and at 0.40 m it is 3,456. That is 1 to 2.5 to 4. Differentiate $p_0 - p_\infty = \rho\Omega^2(r^2 - a^2)$ and the slope left is $2\rho\Omega^2 r$. With $\rho = 1.2$ and $\Omega = 60$, $2\rho\Omega^2$ is 8,640 Pa/m². Times 0.10, 0.25, and 0.40 it reprints 864, 2,160, and 3,456. The central difference landed on the analytic slope.

The 1,080 Pa hole uses the same product. $\rho\Omega^2 a^2 = 1.2 \times 3600 \times 0.25$. Edge dynamic pressure at 30 m/s is $\tfrac{1}{2}\rho u_{\max}^2 = 540$ Pa, half the center deficit. At $r = 0$ the speed is zero, so the static deficit and the total deficit are the same number. The 0.25 m row is 100,380 Pa static and 100,515 Pa total. The 135 Pa gap is the dynamic pressure of 15 m/s. Crossing from 0.10 m to 0.40 m costs $\rho\Omega^2(0.16 - 0.01) = 648$ Pa. Along one radius the gap is zero, because that path is a streamline.

Outer dissipation dies as $1/r^4$. The table gives 0.0313 W/m³ at 0.60 m and 0.0040 W/m³ at 1.00 m, a ratio near 7.8. Shear in the free vortex goes as $1/r^2$, and $\phi = \mu(2S_{r\theta})^2$ squares it. $(1/0.60)^4$ is about 7.7. The core entry 0.0000 is not that curve running inward. Strain is identically zero there.

Dissipation was nonzero where the total pressure was flat#

The last column is flipped. Dissipation is zero inside the core, where the total pressure drops 1,080 Pa, and nonzero outside, where the total pressure is exactly flat.

The reason is that dissipation is attached to strain, not to rotation. Solid-body rotation only turns a fluid element; it does not distort it. The strain-rate tensor is zero, so $\phi = 0$. The free vortex is the opposite. Its vorticity is zero, but $u_\theta \propto 1/r$ means the inner and outer edges sweep past at different speeds. The element is sheared without pause.

The square in the simulation above shows this directly. Inside the core (red) the square keeps its shape as it turns, and outside (purple) it collapses into a parallelogram. What is conserved in a rotating frame and what is not came up in Coriolis force and rothalpy, and the axis is the same here. Spinning and doing work are entered in different ledgers.

Five flows on one table#

Vorticity and dissipation are independent. All four combinations are real.

What matters is that only one cell in the table is bold. The one thing that can be called a loss is $p_0$ decreasing along a streamline. In the other four rows, the spatial variation of $p_0$ is entirely the geometry of the vorticity field, not energy that went missing.

When total-pressure contours look alike for two different reasons#

In practice this distinction breaks at the total-pressure contour on the outlet plane. Feed in an atmospheric boundary layer or a developed pipe profile as the inlet condition and $p_0$ varies in $y$ from the first cell onward. Loss is still zero. The total-pressure deficit a wake creates draws the same contour. That one is real loss.

The two channels already run above, with match outlet spread forcing the outlet $p_0$ spans to match. The outlet contour alone cannot tell them apart. What still separates them is whether a particle's color changes while it travels. The top channel keeps its color. Only the bottom channel changes color as it passes the body.

So when you compute a loss coefficient with a section-averaged $p_0$ as the reference, a shear inflow returns a nonzero value even though there is no loss. The reference has to be that streamline's inlet $p_0$. Even with a mass-flow-weighted average, inlet and outlet have to be averaged the same way for the difference to be loss alone.

One more thing shows up numerically. On a coarse grid, $p_0$ goes artificially flat in rotational regions. When numerical diffusion smears the vorticity, the right side of $\nabla p_0 = \rho\,\mathbf{u}\times\boldsymbol{\omega}$ shrinks, and the core deficit comes out shallower than it really is. On a grid that cuts through a vortex core, the depth of the total-pressure deficit works as an indicator of vorticity resolution. Cut the forced core and the free exterior of one Rankine vortex on a coarse mesh and the 1,080 Pa deficit is the first number to go shallow.

d'Alembert ran into the same place in 1752#

Johann never came around to Newton's theory of viscosity, and neither did his son Daniel nor his student Euler. Solve the flow around a body with inviscid theory alone and the drag comes out zero. That is why the field panicked when d'Alembert published this in 1752. The problem d'Alembert hit with the wave equation in the same period shares one root with it. How wide a class of solutions do you let the equation admit?

Restated in the language we use today, it reads like this. In irrotational inviscid flow $p_0$ is a single constant over the whole domain, so the pressure fore and aft of the body is symmetric and the integral comes to zero. To produce drag, $p_0$ has to fall along a streamline somewhere. What makes that place is viscosity and the vorticity it produces at the wall.

So the question to ask when a total-pressure contour is open in front of you is not "how far did it drop." It is "did it drop along a streamline, or across one?" The first is loss, the second is vorticity.