---
jupyter:
  description: 'The Venturi effect with planex2d and NGSolve: a planar Venturi channel with a 21° convergent and a 7°
    divergent, Navier-Stokes flow at rising Reynolds numbers, Bernoulli''s equation checked along the centerline, and
    the discharge coefficient of the meter from its wall taps.'
---

# The Venturi meter

Why does the pressure drop where a pipe narrows? Where the flow speeds up, part of its pressure has become kinetic
energy, and Bernoulli's equation, $p + \tfrac{1}{2}\rho u^2 = \text{const}$ along a streamline, says how much. This is
the Venturi effect, and the Venturi meter turns it into a flow measurement. From the pressure difference $\Delta p$
between the pipe and the throat, continuity and Bernoulli's equation give the flow rate

$$Q_{ideal} = A_2 \sqrt{\frac{2\,\Delta p}{\rho\,\left(1 - (A_2/A_1)^2\right)}}.$$

Bernoulli's equation neglects viscosity, so the real flow is somewhat less, by the discharge coefficient
$C = Q / Q_{ideal}$. Here, netgen meshes a planar Venturi channel and NGSolve solves the Navier-Stokes equations in it
at rising Reynolds numbers. Bernoulli's equation is checked along the centerline, and the meter's discharge coefficient
comes from its wall taps.

![The speed in a planar Venturi channel at Re 800 with its contour lines, fastest in the throat](../_static/gallery/venturi-meter.png){.result}

```python
import math

import matplotlib.pyplot as plt
import ngsolve as ngs
import numpy as np
import planex as px
from ngsolve.solvers import Newton
```

```python nbsphinx="hidden"
from docs_figures import MAP, Chart, Row


def sampled(cf, label):
    """a field of the flow for Row: cf at points inside the channel"""
    def field(face, plot=True, **_):
        (u0, v0), (u1, v1) = face.bounding_box
        gu, gv = np.meshgrid(np.linspace(u0, u1, 601), np.linspace(v0, v1, 81))
        uv = np.c_[gu.ravel(), gv.ravel()]
        uv = uv[face.Distance(uv, signed=True) < -0.01]
        return {"points": [tuple(p) for p in uv], "values": np.asarray(cf(mesh(uv[:, 0], uv[:, 1]))).ravel(),
                "label": label}
    return field
```

## The channel

The channel follows the classical Venturi tube: a convergent of 21°, a throat as long as it is wide, and a slowly
widening divergent of 7°, which turns part of the speed back into pressure. It is planar, a width of $D = 2$ narrowing
to $d = 1$, so the area ratio is $A_2/A_1 = d/D = 0.5$, where a round tube would have the square of the diameter
ratio. Straight pieces of 4 before and 6 after keep the inlet and the outlet away from the contraction. The upper wall
is a polyline, the lower one its mirror image, and the eight corners between the pieces are rounded with a radius of
1, as the transitions of a real tube are (Figure 1):

```python
conv = 0.5 / math.tan(math.radians(21 / 2))
div = 0.5 / math.tan(math.radians(7 / 2))
x = np.cumsum([0, 4, conv, 1, div, 6])
upper = [(x[0], 1), (x[1], 1), (x[2], 0.5), (x[3], 0.5), (x[4], 1), (x[5], 1)]
channel = px.Face([px.FromPolyline(upper + [(h, -v) for h, v in reversed(upper)], closed=True)])
channel = channel.Filleted([(h, s * v) for h, v in upper[1:5] for s in (1, -1)], 1)
channel.edges.name = "wall"
channel.edges(near=(x[0], 0)).name = "inlet"
channel.edges(near=(x[5], 0)).name = "outlet"
```

```python nbsphinx="hide-input"
Row((channel, "D = 2 narrowing to d = 1, convergent 21°, divergent 7°"), caption="The planar Venturi channel")
```

## The flow

The inlet carries a parabolic profile with a mean velocity of 1, the walls hold the fluid at rest, and the outlet is
left free, so the pressure there is zero on average. With the density 1, the Reynolds number of the throat is
$Re = u_2 d / \nu = 2 / \nu$. Newton's method solves the steady Navier-Stokes equations, each Reynolds number starting
from the solution of the one before:

```python
mesh = ngs.Mesh(px.OCCGeometry(channel).GenerateMesh(maxh=0.08))
X = ngs.VectorH1(mesh, order=2, dirichlet="wall|inlet") * ngs.H1(mesh, order=1)
(u, p), (v, q) = X.TnT()
flow = ngs.GridFunction(X)
flow.components[0].Set(ngs.CF((1.5 * (1 - ngs.y**2), 0)), definedon=mesh.Boundaries("inlet"))


def solve(Re):
    nu = 2 / Re
    form = nu * ngs.InnerProduct(ngs.grad(u), ngs.grad(v)) + ngs.InnerProduct(ngs.grad(u) * u, v) \
        - ngs.div(v) * p - ngs.div(u) * q
    Newton(ngs.BilinearForm(form * ngs.dx), flow, freedofs=X.FreeDofs(), maxerr=1e-10, inverse="umfpack",
           printing=False)


print(mesh.ne, "elements")
```

## Bernoulli along the centerline

The pressure taps sit 1 before the convergent and in the middle of the throat. Along the centerline, the total
pressure $p + \tfrac{1}{2}u^2$ would stay the same without viscosity. Its loss from the inlet tap to the throat,
measured against the throat's dynamic pressure, shrinks as the Reynolds number grows. At the walls, the pressure
difference gives the meter's flow rate $Q_{ideal}$, and the true flow rate is 2, the inlet's width times its mean
velocity:

```python
velocity, pressure = flow.components
inlet_tap, throat_tap = x[1] - 1, (x[2] + x[3]) / 2


def total(s):
    return pressure(mesh(s, 0.0)) + velocity(mesh(s, 0.0))[0] ** 2 / 2


for Re in (20, 100, 200, 400, 800):
    solve(Re)
    loss = (total(inlet_tap) - total(throat_tap)) / (velocity(mesh(throat_tap, 0.0))[0] ** 2 / 2)
    dp = pressure(mesh(inlet_tap, 1 - 1e-6)) - pressure(mesh(throat_tap, 0.5 - 1e-6))
    q_ideal = 1 * math.sqrt(2 * dp / (1 - 0.5**2))
    print("Re", Re, "| Bernoulli loss", round(100 * loss, 1), "% | wall Δp", round(dp, 3), "| C", round(2 / q_ideal, 3))
```

The loss about halves with each doubling of the Reynolds number, so the centerline approaches Bernoulli's equation as
viscosity matters less and less in the core of the flow. The wall taps still read more than Bernoulli's equation gives for the mean
velocities, $\tfrac{1}{2}(2^2 - 1^2) = 1.5$. Near the walls, the fluid is slowed, so the core runs faster than the
mean, and the pressure falls further than the mean velocity says. The meter therefore reads too much flow, and its
discharge coefficient rises toward 1 as the layers near the walls grow thinner. Real meters run at Reynolds numbers far
higher than this laminar flow and are calibrated there.

## The pressure along the channel

At $Re = 800$, the pressure on the centerline falls through the convergent, is lowest in the throat, and rises again in
the divergent (Figure 2). In this laminar flow, the divergent recovers only part of the drop, since the friction at the
walls takes much of it:

```python
recovered = (pressure(mesh(x[4] + 0.2, 1 - 1e-6)) - pressure(mesh(throat_tap, 0.5 - 1e-6))) / dp
print("the divergent recovers", round(100 * recovered), "% of the drop")
```

```python nbsphinx="hide-input"
along = np.linspace(x[0] + 0.1, x[5] - 0.1, 400)
Chart("x", "pressure", curves={"on the centerline": (along, [pressure(mesh(s, 0.0)) for s in along])},
      phases=[(x[1], "convergent"), (x[2], "throat"), (x[3], "divergent")],
      title="The pressure along the channel at Re 800")
```

The speed is highest in the throat, and the layers of slower fluid near the walls grow along the divergent (Figure 3):

```python nbsphinx="hide-input"
Row((channel, "The speed at Re 800", {"field": sampled(ngs.Norm(velocity), "speed"), "cmap": MAP.reversed(),
                                       "clim": (0, 2.5), "bar": "right"}))
```

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from gallery import save_figure
fig, ax = plt.subplots(figsize=(6.4, 2.6))
px.PlotField(velocity, mesh, channel, ax=ax, colormap=MAP.reversed(), min=0, max=2.5, samples=900,
             contours=[0.5, 1.0, 1.5, 2.0])
ax.set_aspect("equal", adjustable="box")
ax.set_xlim(x[1] - 1.5, x[3] + 4)
ax.set_ylim(-1.1, 1.1)
save_figure(fig, "venturi-meter", bare=True, bars={1: "speed"})
```