一つの積分値が0.1406と56.91に同時に読めた — ねじり応力関数とダクト層流

断面上でラプラシアンが-1となる問題を一度解けば、その積分値がねじり定数であり、同時にダクトのf·Reです。構造ソルバと流体ソルバは同じ行列を二度組んでいます。

角棒をねじる問題と角ダクトに水を流す問題

角断面の鋼棒を1 mあたり0.01 radねじるには、トルクがどれだけ要るでしょうか。同じ形のダクトに水を流したとき、摩擦係数はいくつになるでしょうか。この二つは別の学科で別の教科書から習います。ところが答えは、同じ一つの積分値から出てきます。

この記事では、その積分値を実際に計算します。線形三角形要素で断面上のポアソン方程式を一度解き、同じ解を二度読みます。一度はねじり定数 $J$ として、もう一度は層流摩擦群 $f\cdot Re$ として読みます。正方形断面で出るべき値はそれぞれ0.1406と56.91です。どちらもハンドブックに載っている数字です。

3次元問題が断面上のスカラー一つに縮む過程#

まずねじりから見ます。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も合っています。同じ数字なので、間違えようがありません。