A flow pinned between two guesses

The cheapest shape the walls allow showed that the parabola in a pipe is not merely what the equations give. Of every profile that could carry that flow between those walls, it is the one that destroys the least energy — Helmholtz’s minimum-dissipation theorem, computed rather than cited. That essay used the theorem to explain a shape. This one uses it to measure something: the flow through a duct whose shape has no formula.

The idea is old and still underused. A minimum principle does not only single out the true solution. It ranks every other candidate. Any candidate is worse than the true one in a known direction, so every guess is a bound. The principle has a twin, stated in terms of stress rather than velocity, whose guesses are wrong in the other direction. With both, a calculation brackets its own answer. It needs no exact solution to compare with, and no experiment, to say how close it is.

Fully developed flow in any duct #

A long straight duct of any cross-section carries a steady laminar flow along its length. Far from the entrance — past the distance how far before a duct forgets computes — the velocity w depends only on position in the cross-section. It satisfies one equation:

$$-\nabla^2 w = \frac{G}{\mu} \quad\text{in the section}, \qquad w = 0 \text{ on the walls},$$

where G is the pressure gradient along the duct and μ the viscosity. Scale G/μ to one and lengths to a, half the side of the square that will be the test case. The flow rate is Q, the integral of w over the section. Multiplying the equation by w and integrating by parts gives a second expression for Q: it equals the integral of $|\nabla w|^2$. That second integral is the dissipation. In a fully developed duct, the flow rate and the rate of energy loss are the same number — the energy balance that the price of a gradient writes out in general.

Three sections have closed forms: the circle, the equilateral triangle, and the ellipse. The rectangle has a series, Boussinesq’s, which for the square gives Q = 0.56230806 in these units. Almost everything else — an L-shaped duct, a ribbed channel, the passage between the tubes of a heat exchanger — has neither. For those an approximate solution is all there is, and the question is how approximate it is.

Guessing low: the velocity #

Take any trial velocity v that vanishes on the walls and scale it so its dissipation matches its flow. The minimum principle then says that its flow rate cannot exceed the true one:

$$Q \;\ge\; \frac{\left(\int v\right)^2}{\int |\nabla v|^2}.$$

The true flow is the cheapest way to carry its own flow rate. Any other profile, adjusted to cost the same, carries less.

The simplest trial for the square is the product $(1 - x^2)(1 - y^2)$, which is zero on all four walls. Its integrals are elementary — 16/9 and 256/45 — and they give Q ≥ 5/9 = 0.5556. That is 1.2 per cent below the truth, from one line of algebra. Multiply by $(1 + c(x^2 + y^2))$ and choose c to make the bound as large as possible. At c = 0.2025 it reaches 0.56157, 0.13 per cent below the truth. Every value of c gives a valid bound, and none crosses the true value.

Guessing high: the stress #

The twin principle works with stresses instead. The shear stress in the cross-section is a vector field τ. In the true flow it is $\nabla w$, and it must balance the pressure gradient: $\nabla\cdot\tau = -1$. Take any vector field that satisfies that balance — it does not need to come from any velocity, and it does not need to satisfy anything on the walls. Then

$$Q \;\le\; \int |\tau|^2 .$$

The proof is two lines. Write $\tau = \nabla w + \sigma$. Then σ has zero divergence, and the cross term vanishes: $\int\sigma\cdot\nabla w$ integrates by parts to $-\int w\,\nabla\cdot\sigma$ plus a boundary term that is zero because w vanishes on the walls. So $\int|\tau|^2 = Q + \int|\sigma|^2$, and the extra is never negative. This is Kelvin’s principle, older than Helmholtz’s, and it is the complementary energy of elasticity under another name.

The trial stress needs no wall condition. That is what makes the twin easy to use. The simplest balanced field is τ = −(x, y)/2, whose divergence is −1 everywhere. On the square it gives Q ≤ 2/3. Mixing the two directions as −(αx, (1 − α)y) gives a bound for every α, best at α = ½, which is the same field. The bracket from two lines of algebra is therefore 0.5616 to 0.6667. It is lopsided: a smooth velocity is a far better guess than a linear stress, and the lower bound is tighter by more than a factor of a hundred. But it is a bracket, and it came from nothing but the equation.

Both bounds on one mesh #

Hand trials run out quickly. The systematic version gives each principle a finite-element mesh and lets it find the best trial there.

For the velocity, the trial space is the continuous piecewise-linear functions on a mesh of triangles, vanishing on the walls. The best trial in that space is the ordinary finite-element solution, and its flow rate is a guaranteed lower bound. For the stress, every balanced field can be written as $\tau_0 + \nabla^\perp\psi$, where $\tau_0 = -(x, y)/2$ is the particular solution and $\nabla^\perp\psi$ is the rotated gradient of any function ψ. With ψ continuous and piecewise linear on the same mesh, $\nabla^\perp\psi$ is divergence-free exactly — its normal component is continuous across every edge, because it is ψ’s derivative along the edge. The best ψ solves a second finite-element problem, with free edges instead of fixed ones, and gives a guaranteed upper bound.

So one mesh gives two numbers, and the answer lies between them. On the square, at 4, 8, 16, 32 and 64 cells per half-side, the series value 0.5623081 is inside every bracket. At 64 cells the bracket is 0.5621965 to 0.5623592, a width of 0.029 per cent. The calculation certifies that width itself. The series is only a check that the certificate is honest.

A duct with no formula #

The L-shaped duct is the square with one quadrant removed: three unit squares in an L, the kind of passage a bent sheet makes against a wall, or the cross-section of an angle iron. Nobody has a formula for its flow. The two bounds on the same meshes close on it: 0.1891 to 0.2314 at four cells per half-side, 0.2118 to 0.2159 at sixteen, and 0.213991 to 0.214154 at 128 — the flow known to 0.08 per cent without anything to compare it against.

In the engineering form, with the hydraulic diameter 4A/P = 1.5 and the Darcy friction factor, the L’s Poiseuille number lies between 63.04 and 63.09. A number that is only the shape of the hole lists 64 for the circle and 56.91 for the square, and puts the error of treating ducts alike through the hydraulic diameter at a third. The L sits near the circle by that measure, although it looks nothing like one, which says more about the hydraulic diameter than about the L.

What a certificate buys that a convergence study does not #

The usual way to trust a numerical answer is to refine the mesh and watch it settle. On the L that practice is quietly misleading. The lower bound alone reads 0.2066, 0.2118, 0.2134, 0.2138 and 0.2140 as the mesh is halved four times. The differences shrink by a factor of about three each time, and anyone extrapolating would assume the square’s rate of four. Richardson’s extrapolation on that assumption gives 0.21404, which understates the remaining gain, because the true rate is slower; and nothing in the sequence itself says so. The upper bound, run alongside, removes the guesswork. The answer is between the two, whatever the rate, and the width is the error.

A designer needs exactly that. A heat exchanger is sized so the pump can overcome its pressure drop at the worst flow it must carry. What matters there is a number the drop is certain not to exceed, not a best estimate that might sit either side. The velocity bound gives a flow rate that is certainly too small for a given pressure gradient. Turned round, that is a pressure drop certainly too large for a given flow, which is the safe side of the design. The stress bound gives the other side, and says how much margin is being paid for. Both come from one mesh, at the cost of solving one more linear system of the same size.

The two principles also say what each guess must respect and what it may ignore. The velocity must stick to the walls, and nothing else is asked of it. It need not satisfy the momentum equation anywhere, only carry flow, which is the constraint mass has nowhere to go states. The stress must balance the pressure gradient everywhere, and nothing else is asked of it: it may slip along the walls freely. Each principle relaxes exactly the condition the other enforces. That is why their errors fall on opposite sides, and why the true flow, which obeys both, is the only one they agree on.

The corner that slows the bracket #

The two brackets do not close at the same rate. The square’s width falls by a factor of four each time the mesh spacing is halved — rates of 1.99, 2.00 and 2.00 — which is what linear elements give on a smooth solution. The L’s falls by 3.2, then 3.1, then 2.9: rates of 1.67, 1.61 and 1.54, and 1.48 on the finest halving. Those rates are falling towards 4/3.

The reason is the L’s inside corner, where two walls meet at 270° measured through the fluid. Near a corner the solution is fixed by the angle alone. Nothing turns a sharp corner gives the general rule: the flow grows from a corner as $r^n$ with the exponent set by π over the angle. For a 270° corner that exponent is 2/3. The velocity rises as $r^{2/3}$, so its gradient — the shear stress on the two walls that meet there — goes as $r^{-1/3}$ and is infinite at the corner itself. Piecewise-linear elements cannot represent that, and the part of the error the corner contributes shrinks only as $h^{4/3}$ instead of $h^2$. On coarse meshes the smooth part of the error dominates and the rate looks nearly two. As the mesh is refined the corner’s share takes over, and the rate falls towards 4/3.

The velocity on the bisector shows the exponent directly. Fitted between 0.02 and 0.1 half-sides from the corner, it is 0.660, against a theoretical 2/3. The small shortfall is the next term in the corner’s expansion, which is about three per cent of the velocity at a tenth of a half-side and bends the line down there.

That is a limit of the model, and the bracket is where it shows. A real duct’s inside corner is rounded, if only slightly, and the stress there is large but finite. The infinite stress is a statement about a sharp corner that no manufactured duct has. The bracket closing slowly is the calculation reporting that its geometry contains one.

Twisting a shaft is the same problem #

The equation $-\nabla^2 w = 1$ with w = 0 on the boundary is not only duct flow. Saint-Venant’s theory of a twisted shaft reduces to exactly the same equation for a stress function across the shaft’s section. The torsional stiffness is four times the duct’s flow rate, and the shear stress in the shaft is the gradient of the stress function, just as the fluid’s is the gradient of w. Prandtl’s membrane analogy of 1903 — a soap film stretched over a hole the shape of the section and inflated slightly — solves both at once, and engineers used it to measure the torsional stiffness of odd sections before computers could.

So the L’s bracket is also the torsional stiffness of a thick angle section, 0.8560 to 0.8566 in units of $a^4$. Its infinite corner stress is the reason angle irons crack at the inside of the bend and are made with a fillet there. Saint-Venant conjectured, and Pólya proved in 1948, that of all sections of a given area the circle is the stiffest in torsion. It therefore also carries the most flow for its area. The ranking in the figure is that theorem, with the L at the bottom: 0.02378 to 0.02380, against the circle’s 1/8π = 0.03979.

The twin principles have the same pedigree. The velocity principle is the potential energy of elasticity, and the stress principle is the complementary energy — Castigliano’s. Engineers have bounded torsional stiffness from both sides since the 1920s, and the methods of Trefftz and of Pólya and Szegő’s isoperimetric inequalities were largely built on the torsion problem.

What the picture cannot show #

The flow is fully developed. Every number here is for a duct long enough that the entrance is forgotten. Near the entrance the velocity depends on distance along the duct too, the problem is not Poisson’s, and the principles in this form do not apply.

It is laminar. The minimum principle is a property of the Stokes equations, and a duct flow that has become turbulent is no longer the cheapest flow available to it. The cost of going turbulent measures how far it sits from that minimum. Laminar flow in a duct of this size and shape survives to a Reynolds number of about two thousand.

The corner is sharp. The infinite stress and the slow closure of the bracket both belong to an idealised corner. A fillet of radius r replaces the singularity with a stress of order $r^{-1/3}$ and restores the square’s rate on meshes finer than the fillet.

And the bounds are bounds on one number. They pin the flow rate, which is also the dissipation, to within the bracket. They say nothing directly about the velocity at a point, which can be in error by much more than the flow rate — near the corner, by a great deal.

The convention the numbers depend on #

Lengths are in units of a, half the side of the square, so the square is $[-1, 1]^2$ and the L is the square without the quadrant x > 0, y < 0. G/μ is one, so flow rates are in units of $G a^4/\mu$ and velocities of $G a^2/\mu$. The mesh is a grid of squares, each cut on the same diagonal, n per half-side. The Poiseuille number uses the Darcy friction factor and the hydraulic diameter, 2 for the square and 1.5 for the L. The torsional stiffness is four times the flow rate, in units of $a^4$.

How each number was checked #

The square’s series value lies inside the bracket at 4, 8, 16 and 32 cells per half-side, and that containment is the check that the two principles have been implemented the right way round. A sign error in either would put the series outside at once. The square’s bracket must close at a rate of two per halving, to within five hundredths, and does. The L’s must close more slowly than the square’s, faster than 4/3, and more slowly with each halving, and does. The velocity exponent at the corner must be within five hundredths of 2/3. The calculation refuses a mesh of no cells, a fractional mesh, a shape it does not know, and an exponent tolerance of zero.

Who found it, and when #

Kelvin stated the stress principle for a related problem in 1849, and Helmholtz the velocity principle for viscous flow in 1868. Rayleigh used both to bound physical constants in The Theory of Sound. Trefftz in 1927 turned the pair into a method for bounding torsional stiffness, the first systematic use of two-sided bounds in engineering. Prager and Synge in 1947 gave it the geometric form — the true solution lies on a sphere determined by any velocity trial and any stress trial — which is the basis of modern error estimation for finite elements. Boussinesq’s series for the rectangle is from 1868 as well, and Pólya’s proof of Saint-Venant’s conjecture from 1948.

Still open: bounds on the corner itself, and on a duct with an obstacle #

The bracket here is on a global number, the flow rate. The next calculation would bound something local. The wall shear stress near the L’s corner can be written as a known singular function times a coefficient — the stress intensity factor, in the language of fracture mechanics — and the same two principles, with the singular function added to both trial spaces, bound that coefficient from both sides. That would restore the square’s rate of convergence on the L and give the corner’s stress to as many figures as the flow.

Beside it is a duct with something in it: flow past a rod held along the axis of a pipe, or through the passages of a tube bundle. There the walls are not all one boundary, the flow is not simply connected, and the stress principle needs one more condition for each hole — a statement about circulation, which the reciprocal theorem and the circulation essays of this collection supply in other settings. Bounding the pressure drop through a heat exchanger’s tube bundle from both sides is the practical form of that question.