同一个积分值被同时读作 0.1406 和 56.91 — 扭转应力函数与管道层流

在截面上求解一次拉普拉斯算子等于 -1 的问题，得到的积分值既是扭转常数，也是管道的 f·Re。结构求解器和流动求解器在重复组装同一个矩阵。

扭一根方棒的问题，和让水流过方形管道的问题

要把方截面钢棒每米扭转 0.01 rad，需要多大扭矩。同样形状的管道里通水，摩擦系数是多少。这两个问题在不同的院系、用不同的教材讲授。可是答案出自同一个积分值。

本文把这个积分值直接算出来。用线性三角形单元在截面上求解一次泊松方程，再把同一个解读两遍。一遍读作扭转常数 $J$ ，一遍读作层流摩擦群 $f\cdot Re$ 。正方形截面上应当出现的值分别是 0.1406 和 56.91。两个都是手册里的数字。

三维问题收缩成截面上一个标量的过程

先看扭转。Prandtl 没有直接去解应力分量。他改为引入应力函数（stress function，求导之后即为应力的标量场） $\Phi$ 。

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

$\tau_{zx}$ 、 $\tau_{zy}$ 是作用在截面上的剪应力的两个分量。这样设定之后，平衡方程自动满足。剩下的只有一个协调条件，而它正好化为截面 $A$ 上的泊松方程。

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

$G$ 是剪切模量， $\theta$ 是单位长度的扭转角。侧面是自由表面，没有剪应力，因此边界上的 $\Phi$ 是常数。对实心截面可以把这个常数取为 0。扭矩由截面积分取回。

T = 2\int_A \Phi \, dA

接下来是流动。等截面管中流动一旦充分发展，速度就只剩轴向分量 $u(x,y)$ 。 $u$ 沿轴向不变，对流项整项消失。剩下的 Navier–Stokes 是线性的。

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

$\mu$ 是动力黏度， $dp/dz$ 是轴向压力梯度，在截面上为常数。流量也是同样形式的积分。

Q = \int_A u \, dA

两个式子只差符号。在下面的模拟里可以亲手调一调。

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.

拖动长宽比滑块，松弛会从头重新跑一遍，右侧两张卡片会用同一个场填入各自的值。两个按钮并不改变计算。它们只换标签。

逐个符号互换的对照表

扭转	管道层流	共同点
应力函数 $\Phi$	轴向速度 $u$	未知标量
$2G\theta$	$(-dp/dz)/\mu$	常数右端项
自由表面 $\Phi = 0$	无滑移条件 $u = 0$	狄利克雷边界
剪应力 $\lvert\nabla\Phi\rvert$	壁面剪切 $\mu\lvert\nabla u\rvert$	边界上的梯度
扭矩 $T = 2\int \Phi\, dA$	流量 $Q = \int u\, dA$	截面积分
扭转常数 $J$	$f\cdot Re$	由截面形状决定的常数

先做一次归一化会方便很多。求解 $\nabla^2 w = -1$ 、边界上 $w = 0$ 的问题，并记 $I = \int_A w \, dA$ 。于是两个常数这样落出来。

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

前一式由 $\Phi = 2G\theta\, w$ 代入 $T = G\theta J$ 得到。后一式是把 Darcy 摩擦系数 $f = 2 D_h(-dp/dz)/(\rho \bar u^2)$ 与 $Re = \rho \bar u D_h/\mu$ 相乘、消去平均速度 $\bar u$ 的结果。 $D_h = 4A/P$ 是水力直径， $P$ 是湿周。

代入圆截面就是一次验算。半径 $R$ 时 $w = (R^2 - r^2)/4$ 、 $I/A = R^2/8$ 、 $D_h = 2R$ 。代入后 $f\cdot Re = 64$ 原样出现。圆管里指数 4 从何而来已另文讨论过。

一个线性三角形造出的 3×3#

弱形式两边是一样的。乘上试验函数做分部积分后，刚度矩阵里只剩形函数梯度的内积。线性三角形的梯度在单元内是常数。所以不需要积分点，矩阵有闭式表达。

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

对节点 $i,j,k$ 有 $\beta_i = Y_j - Y_k$ 、 $\delta_i = X_k - X_j$ ，其余由下标轮换得到。 $A_e$ 是三角形面积。载荷项按面积的 1/3 均分到三个节点，原因是右端项为常数。用伽辽金加权余量到达同一个矩阵的路径此前已经整理过。

Python — 一次 CG，两个常数#

下面的代码把矩形截面切成结构化网格，每个单元里放两个三角形。用对角预处理 CG 求解一次，再把得到的解读两遍。级数解用于验证。

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

[C] 的最后一行就是本文的要点。50 mm 钢棒的 $J/a^4$ 和 5 mm 水管道的 $4\mu Q/(G_p a^4)$ 是小数点后八位都相同的数。一个给出 702.3 N·m，另一个给出每秒 4.39 mL。

单元尺寸减半时误差移动的比例

[A] 的误差依次是 4.898 → 1.256 → 0.317 → 0.079 %。网格每减半一次，误差约缩小 4 倍。线性三角形的解是 $O(h^2)$ ，它的积分也跟着同样的阶。

符号更有意思。四套网格都把 $J$ 看低、把 $f\cdot Re$ 看高。这不是偶然。有限元位移解总是比精确解更刚。刚度被高估，同样扭矩下就扭得更少，积分值随之变小。换成流动的说法，就是流量被低估。流量偏小，摩擦系数就偏大。同一个偏差在一边读作偏安全，在另一边也读作偏安全，这种情形并不多见。

CG 迭代数从 9 → 32 → 70 → 142 增长。它几乎与网格一边的节点数成正比。这是条件数按 $h^{-2}$ 增长的泊松问题的典型表现，也正是截面变大之后需要多重网格的地方。

把长宽比推上去，两者分别趋向 0.333 和 96#

[B] 中长宽比 1、2、4、8 的 $J/(wh^3)$ 是 0.1405、0.2285、0.2805、0.3065。材料力学教材的 $\beta_2$ 表给出 0.1406、0.229、0.281、0.307，所以到小数点后第三位都对得上。同一行里 $f\cdot Re$ 是 56.95、62.25、73.02、82.50。级数解是 56.91、62.19、72.93、82.34。

两个常数并肩增大，但极限并不相同。截面变薄时 $J/(wh^3)$ 趋向 $1/3$ ，而 $f\cdot Re$ 趋向平行平板的 96。长宽比为 8 时已经到了 0.3065 和 82.5。

最后一列是最大速度与平均速度之比。正方形上得到 2.0975，手册值是 2.096。越薄越向平行平板的 1.5 下降。长宽比 8 时是 1.63。这一列在扭转一侧没有对应物。因为结构关心的是最大应力，流动关心的是最大速度。

皂膜指出最大应力所在的位置

Prandtl 还留下了不做计算、而用实验读出这个方程的办法。在与截面同形的孔上蒙一层皂膜并轻轻吹气，膜的挠度就是 $\Phi$ ，膜的坡度就是剪应力，膜排开的体积就是扭矩。因为承受均匀压力的膜，其控制方程正是同一个泊松方程。

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.

一边增大长宽比，一边跟着红点看。最大坡度总是落在长边的正中间，靠近角点处则衰减到 0。在角点上，膜同时被两条边拉住，几乎是平的。

这一行结论在工程实务里相当有用。受扭的方棒上，裂纹从长边中央开始，而不是从角点开始。管道里壁面剪切在长边中央最大，在角落几乎为 0。方形管道角落里积沉淀物，角落的腐蚀产物不易被冲走，画的都是同一张图。截面越薄，壁面剪切越平，最大值与平均值之比就越接近 1。

这个对应关系断裂的地方

先在空心截面上断裂。边界多于一条时， $\Phi$ 在每条边界上取不同的常数，还要附加确定这些常数的条件。流动一侧没有这种条件。它不过是多了一道墙的狄利克雷问题。

流动这一侧则是假设先崩。入口段里 $u$ 沿轴向变化，对流项复活； $Re$ 超过 2000 之后，“ $f\cdot Re$ 是常数”这句话本身就不成立。黏度被温度牵着走，或者掺进自由面与浮力，右端项就不再是截面上的常数。

在扭转一侧扮演同样角色的是塑性和翘曲约束。端部被约束的薄壁开口截面，已经走到 St. Venant 假设之外。

留下的条件只有两行。右端项在截面上为常数。边界全部是狄利克雷型。只要这两条还成立，两个问题就是同一个问题。

一直在把同一个矩阵组装两遍

结构代码里的扭转模块，和流动代码里的充分发展流动模块，是同一套组装例程被实现了两次。改变的只有一个右端常数，以及解完之后贴在积分值上的标签。

因此验证也一次做完。手上若有新写的扭转求解器，把正方形截面喂进去，顺手把 $f\cdot Re$ 一起取出来即可。只要出现 56.91，结构一侧的 0.1406 也就是对的。同一个数字，错不了。