적분값 하나가 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도 맞다. 같은 숫자니까 틀릴 수가 없다.