---
jupyter:
  description: 'The hydraulic diameter with planex2d and NGSolve: laminar flow through ducts of five sections and the
    ports of an extruded flat tube, the friction factor times the Reynolds number against the exact solutions, and how
    far each duct is from a pipe of the same hydraulic diameter.'
---

# The hydraulic diameter

Does a duct carry flow like a pipe of the same hydraulic diameter? The hydraulic diameter $D_h = 4A/P$, four times the
section's area over its wetted perimeter, is the diameter for a round pipe and the length that pipe formulas take for
ducts of any other shape. For laminar flow, the friction factor of a pipe is $f = 64/Re$, so $f\,Re = 64$ is what the
hydraulic diameter promises. The exact solutions of laminar duct flow say how far each shape keeps that promise. They
are Saint-Venant's for the equilateral triangle, Boussinesq's series for the rectangle and the closed form for the
annulus.

Here, NGSolve solves the laminar flow through five sections and through the ports of a flat tube, and $f\,Re$ of each
is checked against its exact value.

![The laminar velocity in five duct sections](../_static/gallery/hydraulic-diameter.png){.result}

```python
import math

import ngsolve as ngs
import numpy as np
import planex as px
```

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


def sampled(field, face):
    """the velocity for Row, at grid points inside the section"""
    def values(f, plot=True, **_):
        (u0, v0), (u1, v1) = f.bounding_box
        step = min(u1 - u0, v1 - v0) / 120
        gu, gv = np.meshgrid(np.arange(u0, u1 + step, step), np.arange(v0, v1 + step, step))
        uv = np.c_[gu.ravel(), gv.ravel()]
        uv = uv[f.Distance(uv, signed=True) < -1e-3]
        return {"points": [tuple(p) for p in uv], "values": np.asarray(field(uv)).ravel(), "label": "velocity"}
    return values
```

## Laminar flow in a duct

In fully developed laminar flow, the velocity $w$ along the duct satisfies $\nabla^2 w = -G/\mu$ in the section, with
$w = 0$ on the walls, for the pressure gradient $G$ and the viscosity $\mu$. With $G/\mu = 1$, the flow rate is
$Q = \int w\,dA$, and the friction factor times the Reynolds number is

$$f\,Re = \frac{2 D_h^2 A}{Q}.$$

`HydraulicDiameter()` measures $D_h$, and `duct` solves for the velocity and returns $f\,Re$ with the velocity field:

```python
def duct(face):
    Dh = face.Measure(px.measurements.HydraulicDiameter())
    mesh = ngs.Mesh(px.OCCGeometry(face).GenerateMesh(maxh=0.1))
    mesh.Curve(3)
    V = ngs.H1(mesh, order=4, dirichlet=".*")
    w, v = V.TnT()
    K = ngs.BilinearForm(ngs.grad(w) * ngs.grad(v) * ngs.dx).Assemble()
    F = ngs.LinearForm(v * ngs.dx).Assemble()
    velocity = ngs.GridFunction(V)
    velocity.vec.data = K.mat.Inverse(V.FreeDofs(), inverse="umfpack") * F.vec
    return 2 * Dh**2 * face.mass / ngs.Integrate(velocity, mesh), lambda uv: velocity(mesh(uv[:, 0], uv[:, 1]))
```

## The exact flows

For a rectangle of half sides $a \le b$, Boussinesq's series gives the flow rate
$Q = \tfrac{4}{3} b a^3 \left(1 - \tfrac{192 a}{\pi^5 b} \sum_{i\ \text{odd}} \tanh(i\pi b / 2a)/i^5\right)$. For an
annulus of radii $r_i < r_o$, it is $Q = \tfrac{\pi}{8}\left(r_o^4 - r_i^4 - (r_o^2 - r_i^2)^2 / \ln(r_o/r_i)\right)$,
and the equilateral triangle has $f\,Re = 160/3$ exactly:

```python
def rectangle(width, height):
    a, b = min(width, height) / 2, max(width, height) / 2
    series = sum(math.tanh(i * math.pi * b / (2 * a)) / i**5 for i in range(1, 200, 2))
    Q = 4 / 3 * b * a**3 * (1 - 192 * a / (math.pi**5 * b) * series)
    A, P = width * height, 2 * (width + height)
    return 2 * (4 * A / P) ** 2 * A / Q


def annulus(ro, ri):
    Q = math.pi / 8 * (ro**4 - ri**4 - (ro**2 - ri**2) ** 2 / math.log(ro / ri))
    return 2 * (2 * (ro - ri)) ** 2 * math.pi * (ro**2 - ri**2) / Q
```

## Five sections

The sections are a circle, a square, a rectangle of 4:1, an equilateral triangle and an annulus whose inner radius is
half the outer one. Against the pipe's 64, the triangle and the square carry more flow than their hydraulic diameter
promises, and the flat rectangle and the annulus less (Figure 1):

```python
sections = {
    "circle": (px.Circle((0, 0), 1).Face(), 64),
    "square": (px.QuickSketch().RectangleC(2, 2).Face(), rectangle(2, 2)),
    "rectangle 4:1": (px.QuickSketch().RectangleC(8, 2).Face(), rectangle(8, 2)),
    "triangle": (px.QuickSketch().Polygon((0, 0), 1, 3).Face(), 160 / 3),
    "annulus": (px.Circle((0, 0), 1).Face() - px.Circle((0, 0), 0.5).Face(), annulus(1, 0.5)),
}
flows = {}
for name, (face, exact) in sections.items():
    fRe, flows[name] = duct(face)
    print(name, "| f Re", round(fRe, 4), "| exact", round(exact, 4), "| against a pipe", f"{fRe / 64 - 1:+.1%}")
```

```python nbsphinx="hide-input"
Row(*[(face, name.capitalize(), {"field": sampled(flows[name], face), "cmap": MAP.reversed()})
      for name, (face, exact) in sections.items()], caption="The laminar velocity in the five sections")
```

The corners of the square and the triangle carry hardly any flow, which the wetted perimeter still counts, and the
narrow gap of the flat rectangle and the annulus holds the flow back more than a pipe's wall would. A pipe formula with
the hydraulic diameter is therefore off by up to half for these laminar ducts, in both directions.

## A flat tube

The condensers and radiators of cars carry their refrigerant or coolant in flat aluminum tubes, extruded with a row
of small ports. A tube 16 × 2 mm with rounded ends and walls of 0.35 mm holds eleven ports of 1.05 × 1.3 mm between
webs of 0.3 mm. `HydraulicDiameter` measures the ports together, four times their area over their wetted perimeter,
which is one port's $D_h$. Every port is a rectangle, so $f\,Re$ is Boussinesq's again (Figure 2):

```python
tube = px.QuickSketch((-7, -1)).Line(14).Arc(1, 180).Line(14).Arc(1, 180).Face()
port_w, port_h = (14.6 - 10 * 0.3) / 11, 2 - 2 * 0.35
ports = px.Glue([px.QuickSketch((-7.3 + k * (port_w + 0.3), -port_h / 2)).Rectangle(port_w, port_h).Face()
                 for k in range(11)])
fRe, flow = duct(ports)
print("D_h", round(ports.Measure(px.measurements.HydraulicDiameter()), 4), "mm | f Re", round(fRe, 4),
      "| exact", round(rectangle(port_w, port_h), 4))
```

```python nbsphinx="hide-input"
view = ((-8.2, 8.2), (-1.2, 1.2))
Row((ports, "The velocity in the ports", {"field": sampled(flow, ports), "cmap": MAP.reversed(), "clim": (0, 0.1),
                                          "fills": [(tube - ports, "#c3cbd3")], "view": view}),
    caption="The laminar flow in an extruded flat tube")
```

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
import matplotlib.pyplot as plt
from gallery import save_figure
places = {"rectangle 4:1": (5.1, 2.8), "circle": (1, 0), "square": (3.8, 0), "triangle": (6.4, 0), "annulus": (9.2, 0)}


def all_sections(region, plot=True, **_):
    """the velocity of each section against its own largest, for PlotField on the sections laid out side by side"""
    points, values = [], []
    for name, (face, exact) in sections.items():
        (u0, v0), (u1, v1) = face.bounding_box
        gu, gv = np.meshgrid(np.arange(u0, u1 + 0.025, 0.025), np.arange(v0, v1 + 0.025, 0.025))
        uv = np.c_[gu.ravel(), gv.ravel()]
        uv = uv[face.Distance(uv, signed=True) < -1e-3]
        w = np.asarray(flows[name](uv)).ravel()
        points += [(u + places[name][0], v + places[name][1]) for u, v in uv]
        values += list(w / w.max())
    return {"points": points, "values": np.array(values), "label": "velocity against the largest"}


fig, ax = plt.subplots(figsize=(6.4, 3.4))
region = px.Glue([face.Translated(places[name]) for name, (face, exact) in sections.items()])
px.PlotField(all_sections, region, ax=ax, colormap=MAP.reversed(), min=0, max=1)
save_figure(fig, "hydraulic-diameter", bare=True, bars={1: "velocity against the largest"})
```