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.

import math

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

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:

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:

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

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%}")

circle | f Re 64.0 | exact 64 | against a pipe -0.0%
square | f Re 56.9083 | exact 56.9083 | against a pipe -11.1%
rectangle 4:1 | f Re 72.9311 | exact 72.9311 | against a pipe +14.0%
triangle | f Re 53.3333 | exact 53.3333 | against a pipe -16.7%
annulus | f Re 95.2505 | exact 95.2502 | against a pipe +48.8%

(a) Circle

(b) Square

(c) Rectangle 4:1

(d) Triangle

(e) Annulus

Figure 1: 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):

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

D_h 1.1645 mm | f Re 57.4392 | exact 57.4392

Figure 2: The laminar flow in an extruded flat tube