One Integral Read as 0.1406 and 56.91 — Torsion Stress Function and Duct Laminar Flow

Solve the problem whose Laplacian is -1 over the cross-section once, and that integral is both the torsion constant and the duct's f·Re. Structural solvers and flow solvers have been assembling the same matrix twice.

Twisting a square bar, and pushing water through a square duct#

How much torque does it take to twist a square steel bar by 0.01 rad per meter? Push water through a duct of the same shape, and what is the friction factor? The two questions are taught in different departments out of different textbooks. Yet both answers come from a single integral.

This post computes that integral directly. One Poisson solve on the cross-section with linear triangular elements, then the same solution read twice. Once as the torsion constant $J$ , once as the laminar friction group $f\cdot Re$ . For a square section the values that should come out are 0.1406 and 56.91. Both are handbook numbers.

How a three-dimensional problem shrinks to one scalar on a cross-section#

Start with torsion. Prandtl did not solve for the stress components directly. Instead he set up a stress function (a scalar field whose derivatives are stresses) $\Phi$ .

\tau_{zx} = \frac{\partial \Phi}{\partial y}, \qquad \tau_{zy} = -\frac{\partial \Phi}{\partial x}

$\tau_{zx}$ and $\tau_{zy}$ are the two components of shear stress acting on the cross-section. Defined this way, the equilibrium equations are satisfied automatically. What is left is one compatibility condition, and it becomes a Poisson equation on the cross-section $A$ .

\nabla^2 \Phi + 2G\theta = 0, \qquad \Phi = 0 \ \text{ on } \ \partial A

$G$ is the shear modulus, $\theta$ the twist angle per unit length. The side surfaces are free, so there is no shear there, and $\Phi$ is constant on the boundary. For a solid section that constant can be set to zero. The torque is recovered by a section integral.

T = 2\int_A \Phi \, dA

Now the flow. In a duct of constant cross-section, fully developed flow leaves only the axial component $u(x,y)$ . Since $u$ does not vary along the axis, the convective term disappears entirely. What remains of Navier–Stokes is linear.

\mu \nabla^2 u = \frac{dp}{dz}, \qquad u = 0 \ \text{ on the wall}

$\mu$ is the dynamic viscosity and $dp/dz$ the axial pressure gradient, constant over the cross-section. The flow rate is an integral of the same shape.

Q = \int_A u \, dA

The two sets of equations differ only in their symbols. Try it in the simulation below.

aspect ratio w/h1.0

J/(w h³) 0.00000f·Re 0.000sweeps 0

Drag the aspect ratio: the solve restarts and both cards refill from the one field. Watch that the two numbers move in opposite directions — J/(w h³) climbs toward 1/3 while f·Re climbs toward 96 — and that neither card ever needs a second matrix.

Move the aspect-ratio slider and the relaxation restarts from scratch, with the two cards on the right filling in from the same field. The two buttons do not change the computation. They only change the labels.

A table that swaps the symbols one at a time#

Torsion	Duct laminar flow	Shared
Stress function $\Phi$	Axial velocity $u$	Unknown scalar
$2G\theta$	$(-dp/dz)/\mu$	Constant right-hand side
Free surface $\Phi = 0$	No-slip $u = 0$	Dirichlet boundary
Shear stress $\lvert\nabla\Phi\rvert$	Wall shear $\mu\lvert\nabla u\rvert$	Slope at the boundary
Torque $T = 2\int \Phi\, dA$	Flow rate $Q = \int u\, dA$	Section integral
Torsion constant $J$	$f\cdot Re$	Constant fixed by the section shape

Normalizing once makes life easier. Solve $\nabla^2 w = -1$ with $w = 0$ on the boundary, and let $I = \int_A w \, dA$ . The two constants then drop out like this.

J = 4I, \qquad f\cdot Re = \frac{2 D_h^2}{I/A}

The first follows from putting $\Phi = 2G\theta\, w$ into $T = G\theta J$ . The second comes from multiplying the Darcy friction factor $f = 2 D_h(-dp/dz)/(\rho \bar u^2)$ by $Re = \rho \bar u D_h/\mu$ , which erases the mean velocity $\bar u$ . $D_h = 4A/P$ is the hydraulic diameter and $P$ the wetted perimeter.

Plug in a circular section and it checks out. For radius $R$ , $w = (R^2 - r^2)/4$ , $I/A = R^2/8$ , $D_h = 2R$ . Substituting gives $f\cdot Re = 64$ directly. Where the fourth power in a round pipe comes from was covered separately.

The 3×3 that a single linear triangle produces#

The weak form is the same on both sides. Multiply by a test function, integrate by parts, and the stiffness matrix keeps only the inner product of shape-function gradients. On a linear triangle the gradient is constant inside the element. So no quadrature points are needed and the matrix comes out in closed form.

K^{(e)}_{mn} = \frac{\beta_m \beta_n + \delta_m \delta_n}{4 A_e}, \qquad F^{(e)}_m = \frac{A_e}{3}

For nodes $i,j,k$ , $\beta_i = Y_j - Y_k$ and $\delta_i = X_k - X_j$ , with the rest obtained by cycling the indices. $A_e$ is the triangle area. The load term spreads one third of the area onto each of the three nodes because the right-hand side is constant. The route to the same matrix through Galerkin weighted residuals was written up earlier.

Python — two constants from one CG solve#

The code below cuts a rectangular section with a structured grid and puts two triangles in each cell. It solves once with diagonally preconditioned CG and reads the result twice. The series solution is there for verification.

import math
 
 
def tri_stiffness(p0, p1, p2):
    """Linear triangle: K = (beta_i beta_j + delta_i delta_j) / (4A)."""
    (x0, y0), (x1, y1), (x2, y2) = p0, p1, p2
    a2 = x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1)
    area = 0.5 * a2
    beta = (y1 - y2, y2 - y0, y0 - y1)
    delta = (x2 - x1, x0 - x2, x1 - x0)
    k = [[(beta[r] * beta[c] + delta[r] * delta[c]) / (2.0 * a2)
          for c in range(3)] for r in range(3)]
    return k, area
 
 
def build_mesh(w, h, nx, ny):
    nodes, idx = [], {}
    for j in range(ny + 1):
        for i in range(nx + 1):
            idx[(i, j)] = len(nodes)
            nodes.append((w * i / nx, h * j / ny))
    tris = []
    for j in range(ny):
        for i in range(nx):
            a, b = idx[(i, j)], idx[(i + 1, j)]
            c, d = idx[(i + 1, j + 1)], idx[(i, j + 1)]
            tris.append((a, b, c))
            tris.append((a, c, d))
    fixed = set()
    for j in range(ny + 1):
        for i in range(nx + 1):
            if i in (0, nx) or j in (0, ny):
                fixed.add(idx[(i, j)])
    return nodes, tris, fixed
 
 
def assemble_poisson(nodes, tris, fixed):
    """-lap(u) = 1 with u = 0 on 'fixed'.  Returns CSR-ish rows and rhs."""
    n = len(nodes)
    rows = [dict() for _ in range(n)]
    rhs = [0.0] * n
    for (a, b, c) in tris:
        k, area = tri_stiffness(nodes[a], nodes[b], nodes[c])
        ids = (a, b, c)
        for r in range(3):
            if ids[r] in fixed:
                continue
            rhs[ids[r]] += area / 3.0
            for c2 in range(3):
                if ids[c2] in fixed:
                    continue
                rows[ids[r]][ids[c2]] = rows[ids[r]].get(ids[c2], 0.0) + k[r][c2]
    for f in fixed:
        rows[f] = {f: 1.0}
        rhs[f] = 0.0
    return rows, rhs
 
 
def cg_solve(rows, rhs, tol=1e-12, itmax=20000):
    n = len(rhs)
    x = [0.0] * n
    r = rhs[:]
    z = [r[i] / rows[i][i] for i in range(n)]
    p = z[:]
    rz = sum(r[i] * z[i] for i in range(n))
    r0 = math.sqrt(sum(v * v for v in r))
    for it in range(itmax):
        ap = [0.0] * n
        for i in range(n):
            s = 0.0
            for j, v in rows[i].items():
                s += v * p[j]
            ap[i] = s
        alpha = rz / sum(p[i] * ap[i] for i in range(n))
        for i in range(n):
            x[i] += alpha * p[i]
            r[i] -= alpha * ap[i]
        rn = math.sqrt(sum(v * v for v in r))
        if rn <= tol * r0:
            return x, it + 1
        z = [r[i] / rows[i][i] for i in range(n)]
        rz2 = sum(r[i] * z[i] for i in range(n))
        beta = rz2 / rz
        rz = rz2
        p = [z[i] + beta * p[i] for i in range(n)]
    return x, itmax
 
 
def section_integral(nodes, tris, u):
    tot = 0.0
    for (a, b, c) in tris:
        _, area = tri_stiffness(nodes[a], nodes[b], nodes[c])
        tot += area * (u[a] + u[b] + u[c]) / 3.0
    return tot
 
 
def series_rect(w, h, nterm=60):
    """Exact integral of the Prandtl/duct solution over a w x h rectangle."""
    s = h if h < w else w
    lg = w if h < w else h
    acc = 0.0
    for m in range(1, 2 * nterm, 2):
        acc += math.tanh(m * math.pi * lg / (2.0 * s)) / m ** 5
    j = (1.0 / 3.0) * lg * s ** 3 * (1.0 - (192.0 / math.pi ** 5) * (s / lg) * acc)
    return j / 4.0          # integral of u  ==  J / 4
 
 
def solve_section(w, h, nx, ny):
    nodes, tris, fixed = build_mesh(w, h, nx, ny)
    rows, rhs = assemble_poisson(nodes, tris, fixed)
    u, its = cg_solve(rows, rhs)
    iu = section_integral(nodes, tris, u)
    area, perim = w * h, 2.0 * (w + h)
    dh = 4.0 * area / perim
    return dict(int_u=iu, jtor=4.0 * iu, umean=iu / area,
                fre=2.0 * dh * dh / (iu / area), umax=max(u), its=its,
                ndof=len(nodes))
 
 
if __name__ == '__main__':
    ex_i = series_rect(1.0, 1.0)
    print("[A] square bar, one Poisson solve -> two constants  (exact J/a^4 = %.6f,"
          " f*Re = %.4f)" % (4 * ex_i, 2.0 / ex_i))
    print("  mesh    nodes    integral u      J/a^4    err(%)      f*Re   err(%)  CG")
    for n in (8, 16, 32, 64):
        r = solve_section(1.0, 1.0, n, n)
        print("%3dx%-3d %7d  %.8f  %.6f  %7.3f  %8.4f  %7.3f %4d"
              % (n, n, r['ndof'], r['int_u'], r['jtor'],
                 100 * (r['jtor'] / (4 * ex_i) - 1), r['fre'],
                 100 * (r['fre'] / (2.0 / ex_i) - 1), r['its']))
 
    print()
    print("[B] aspect-ratio sweep, 64x64 mesh   (w x h = AR x 1)")
    print("  AR   J/(w h^3)   beta2(ref)     f*Re   fRe(exact)   u_max/u_mean")
    REF_B2 = {1: 0.1406, 2: 0.2290, 4: 0.2810, 8: 0.3070}
    for ar in (1, 2, 4, 8):
        w, h = float(ar), 1.0
        r = solve_section(w, h, 64, 64)
        ex = series_rect(w, h)
        dh = 4.0 * w * h / (2.0 * (w + h))
        print("%4d   %.5f      %.4f    %8.4f   %8.4f   %8.4f"
              % (ar, r['jtor'] / (w * h ** 3), REF_B2[ar], r['fre'],
                 2 * dh * dh / (ex / (w * h)), r['umax'] / r['umean']))
 
    print()
    print("[C] the same integral read twice (square section, 64x64)")
    r = solve_section(1.0, 1.0, 64, 64)
    print("  dimensionless integral of u over A = %.8f a^4" % r['int_u'])
 
    G, THETA, SIDE = 80e9, 0.01, 0.05          # steel bar, 50 mm square
    j = r['jtor'] * SIDE ** 4
    print("  steel bar  a=50 mm, G=80 GPa, twist=0.01 rad/m")
    print("     J = 4*int*a^4 = %.4e m^4    T = G*theta*J = %.1f N.m" % (j, G * THETA * j))
 
    MU, RHO, DPDX, HALF = 1.0e-3, 1000.0, 200.0, 0.005   # water, 5 mm square duct
    q = DPDX / MU * r['int_u'] * HALF ** 4
    area = HALF ** 2
    ubar = q / area
    dh = HALF
    re = RHO * ubar * dh / MU
    print("  water duct a=5 mm, dp/dx=200 Pa/m, mu=1e-3 Pa.s")
    print("     Q = (G_p/mu)*int*a^4 = %.3e m^3/s   u_mean = %.4f m/s   Re = %.0f"
          % (q, ubar, re))
    print("     f = (f*Re)/Re = %.4f      f*Re = %.4f  (handbook 56.91)"
          % (r['fre'] / re, r['fre']))
    print("  J/a^4 (bar) = %.8f   vs   4*mu*Q/(G_p*a^4) (duct) = %.8f"
          % (j / SIDE ** 4, 4 * MU * q / DPDX / HALF ** 4))

[A] square bar, one Poisson solve -> two constants  (exact J/a^4 = 0.140577, f*Re = 56.9083)
  mesh    nodes    integral u      J/a^4    err(%)      f*Re   err(%)  CG
  8x8        81  0.03342303  0.133692   -4.898   59.8390    5.150    9
 16x16      289  0.03470275  0.138811   -1.256   57.6323    1.272   32
 32x32     1089  0.03503302  0.140132   -0.317   57.0890    0.318   70
 64x64     4225  0.03511638  0.140466   -0.079   56.9535    0.079  142
 
[B] aspect-ratio sweep, 64x64 mesh   (w x h = AR x 1)
  AR   J/(w h^3)   beta2(ref)     f*Re   fRe(exact)   u_max/u_mean
   1   0.14047      0.1406     56.9535    56.9083     2.0975
   2   0.22848      0.2290     62.2469    62.1922     1.9932
   4   0.28047      0.2810     73.0194    72.9311     1.7758
   8   0.30647      0.3070     82.4998    82.3386     1.6315
 
[C] the same integral read twice (square section, 64x64)
  dimensionless integral of u over A = 0.03511638 a^4
  steel bar  a=50 mm, G=80 GPa, twist=0.01 rad/m
     J = 4*int*a^4 = 8.7791e-07 m^4    T = G*theta*J = 702.3 N.m
  water duct a=5 mm, dp/dx=200 Pa/m, mu=1e-3 Pa.s
     Q = (G_p/mu)*int*a^4 = 4.390e-06 m^3/s   u_mean = 0.1756 m/s   Re = 878
     f = (f*Re)/Re = 0.0649      f*Re = 56.9535  (handbook 56.91)
  J/a^4 (bar) = 0.14046553   vs   4*mu*Q/(G_p*a^4) (duct) = 0.14046553

The last line of [C] is the point of this post. The $J/a^4$ of a 50 mm steel bar and the $4\mu Q/(G_p a^4)$ of a 5 mm water duct are the same number to eight decimal places. One yields 702.3 N·m, the other 4.39 mL per second.

How the error moves when the elements are halved#

The error in [A] runs 4.898 → 1.256 → 0.317 → 0.079 %. Each halving of the grid spacing cuts it by about a factor of four. The linear-triangle solution is $O(h^2)$ and its integral follows the same order.

The sign is more interesting. All four grids read $J$ low and $f\cdot Re$ high. That is no accident. A finite element displacement solution is always stiffer than the exact one. Overestimated stiffness means less twist for the same torque, and a smaller integral. Translated into flow language, the flow rate is underestimated. A smaller flow rate makes the friction factor larger. This is one of the rare cases where the same bias reads as conservative on one side and conservative on the other as well.

The CG iteration count grew 9 → 32 → 70 → 142. That is nearly proportional to the node count along one side of the grid. It is the classic Poisson behavior of a condition number growing like $h^{-2}$ , and the point where a larger section starts to need multigrid.

Push the aspect ratio and the two split toward 0.333 and 96#

In [B], $J/(wh^3)$ for aspect ratios 1, 2, 4, 8 comes out as 0.1405, 0.2285, 0.2805, 0.3065. The $\beta_2$ table in a mechanics-of-materials textbook reads 0.1406, 0.229, 0.281, 0.307, so they agree to three decimal places. On the same rows $f\cdot Re$ is 56.95, 62.25, 73.02, 82.50. The series solution gives 56.91, 62.19, 72.93, 82.34.

The two constants grow side by side, but their limits differ. As the section gets thin, $J/(wh^3)$ goes to $1/3$ and $f\cdot Re$ goes to 96, the parallel-plate value. At aspect ratio 8 it has already reached 0.3065 and 82.5.

The last column is the ratio of maximum to mean velocity. The square gives 2.0975, and the handbook value is 2.096. The thinner the section, the closer it drops to the parallel-plate value of 1.5. At aspect ratio 8 it is 1.63. This column has no counterpart on the torsion side. Structures ask about maximum stress, flows about maximum velocity.

A soap film points at where the maximum stress sits#

Prandtl also left a way to read this equation by experiment instead of computation. Stretch a soap film over a hole shaped like the cross-section and blow gently: the deflection of the film is $\Phi$ , the slope of the film is the shear stress, and the volume the film pushes out is the torque. The governing equation of a membrane under uniform pressure is the very same Poisson equation.

aspect ratio w/h2.0

peak on the longc1 0.0000τmax/τmean 0.000

Stretch the section and watch the red marker: it stays at the middle of the long side and never moves to a corner, where the film is flat. Switch the two buttons — the bars are not redrawn, only relabelled.

Stretch the aspect ratio and follow the red dot. The maximum slope always sits at the midpoint of the long side and dies to zero near the corners. At a corner the film is held by two edges at once and stays nearly flat.

In practice that one line is quite useful. In a square bar under torsion, cracks start at the middle of the long side, not at the corners. In a duct the wall shear is maximum at the middle of the long side and nearly zero in the corners. Sediment settling in the corners of a rectangular duct, and corrosion products in the corners not washing away, are the same picture. The thinner the section, the flatter the wall shear, and the closer the max-to-mean ratio gets to 1.

Where this correspondence breaks#

It breaks first on hollow sections. With more than one boundary, $\Phi$ takes a different constant on each boundary, and an extra condition comes along to fix those constants. The flow side has no such condition. It is just a Dirichlet problem with one more wall.

On the flow side the assumptions collapse first. In the entrance region $u$ varies along the axis and the convective term comes back to life, and once $Re$ passes 2000 the statement that $f\cdot Re$ is a constant stops holding at all. If viscosity is dragged around by temperature, or a free surface or buoyancy gets involved, the right-hand side is no longer constant over the cross-section.

The same role on the torsion side is played by plasticity and warping restraint. A thin open section with restrained ends steps outside the St. Venant assumptions.

Two conditions are left. The right-hand side must be constant over the cross-section. Every boundary must be Dirichlet. As long as both hold, the two problems are one problem.

The same matrix, assembled twice#

The torsion module in a structural code and the fully developed flow module in a flow code are two implementations of the same assembly routine. What changes is one right-hand-side constant, and the label attached to the integral once the solve is done.

So verification finishes in one pass too. If you have just written a torsion solver, run it on a square section and pull out $f\cdot Re$ along with the rest. If 56.91 comes out, the 0.1406 on the structural side is right too. It is the same number, so it cannot be wrong.