지름을 1% 잘못 재면 유량은 4% 틀린다 — 벽 반 칸과 Poiseuille의 상수

점성 스킴이 완벽해도 벽을 반 칸 안쪽에 놓으면 유량은 그 네 배로 어긋난다.

검증 케이스가 4% 어긋났는데 스킴은 멀쩡했다#

원관 층류는 해석해가 있는 몇 안 되는 검증 케이스다. 새로 짠 점성항 이산화를 여기에 걸었더니 유량이 4% 낮게 나왔다. 스킴을 의심하는 것이 순서였다. 그런데 격자를 두 배로 조밀하게 해도 오차는 절반으로 줄지 않고 그대로 남았다. 수렴하지 않는 오차는 이산화 오차가 아니라 기하 오차다.

범인은 벽이 서 있는 위치였다. 반지름을 1% 안쪽으로 잘못 잡으면 유량은 4% 줄어든다. 이 글은 그 배율 4가 어디서 오는지, 그리고 같은 지수 4를 1838년 Poiseuille가 유리관 실험에서 어떻게 거꾸로 이용했는지를 따라간다.

유량은 지름의 네 제곱에 걸려 있다#

완전발달 원관 유동에서 남는 것은 축방향 운동량 균형 하나다.

\frac{\mu}{r}\frac{d}{dr}\left(r\,\frac{du}{dr}\right) = \frac{dp}{dx}

$\mu$ 는 점성계수, $r$ 은 중심축에서 잰 반지름 좌표, $u$ 는 축방향 속도다. 중심에서 $du/dr = 0$ , 벽에서 $u = 0$ 을 걸면 포물선 하나가 나온다.

u(r) = \frac{\Delta p}{4\mu L}\left(R^{2} - r^{2}\right)

$\Delta p$ 는 길이 $L$ 구간의 압력 강하, $R$ 은 관 반지름이다. 이것을 단면에 적분하면 유량이 된다.

Q = \int_{0}^{R} u(r)\,2\pi r\,dr = \frac{\pi\,\Delta p\,D^{4}}{128\,\mu L}

$D = 2R$ 은 관 지름이다. 적분에서 $r\,dr$ 이 한 번 더 곱해지면서 $R^{2}$ 가 $R^{4}$ 가 된다. 지름의 네 제곱은 여기서 나온다.

이 지수가 오차에 걸리면 배율이 된다.

\frac{\delta Q}{Q} = 4\,\frac{\delta D}{D} + \mathcal{O}\!\left(\left(\frac{\delta D}{D}\right)^{2}\right)

$\delta D$ 는 지름 오차, $\delta Q$ 는 그로 인한 유량 오차다. 벽 위치 오차는 다른 검증량에서는 보통 1차로 나타난다. 유량에서만 네 배로 갚는다.

아래 시뮬레이션에서 직접 조작해보자.

dD/D-2.5%dp/dx1.0

Drag dD/D to 1 % and watch the red bar: the wall moved by 1 %, the flow rate moved by 4.06 %. The tracer counters at the right are an independent measurement — they know nothing about the formula, and they still settle on 0.904 (currently 1.000 after 0 particles). Switch to wall pinned at cell centre and lower N: a mesh of 10 cells loses 19 % of the flow with a perfectly correct viscous scheme.

dD/D 슬라이더로 아래쪽 관의 벽만 옮기고, 오른쪽 붉은 막대와 카운터를 본다. 카운터는 공식을 모른 채 출구를 지나간 입자를 세는 독립적인 측정인데도 $(1+e)^{4}$ 로 수렴한다. wall pinned at cell centre로 바꾸고 N을 10까지 내리면 점성 스킴이 완벽한 상태에서 유량이 19% 사라진다.

1822년 Navier가 남긴 자리, 1838년 Poiseuille가 채운 지수#

Navier는 교량 엔지니어였다. 그는 유체 저항이 실측으로는 분명한데 오일러 방정식에는 표현이 없다는 점에 붙들려 있었다. 1822년, 아직 탄성학의 수학적 기반이 정리되기 전에 그는 속도의 라플라시안 형태로 점성항을 유도해 오일러 방정식에 붙였다. 달랑베르 이후 반세기 넘게 정체돼 있던 유체역학이 다시 움직인 지점이다.

같은 시기 École Polytechnique에는 Cauchy가 교수로, Coriolis가 강사로 있었다. 회전 좌표계 에너지 방정식이 수차 실험에서 나온 것도 이 학교의 같은 공기 속에서다. 1816년 학교가 잠시 폐쇄되자 입학생 Poiseuille는 진로를 바꿔 의대로 갔다.

그는 의사가 되어 혈류를 쟀다. 혈압이 혈관 지름에 따라 어떻게 달라지는지가 그의 질문이었다. 미세혈관을 흉내내려고 유리관을 직접 뽑았고, 가장 가는 것은 지름 0.015 mm였다. 머리카락 굵기의 5분의 1이다. 1838년 그가 발표한 것은 이런 형태였다.

Q = K''\,\frac{\Delta p\,D^{4}}{L}, \qquad K'' = \frac{\pi}{128\,\mu}

오른쪽의 $K''$ 는 나중에 채워진 값이다. Poiseuille 자신은 점도라는 개념을 갖고 있지 않았다. 그는 $K''$ 를 그냥 상수로 적었다. 실험이 확정한 것은 상수가 아니라 지수였다.

Python으로 재본 "벽 반 칸"의 값#

위 운동량 균형을 반지름 방향 유한체적으로 풀었다. 벽 무활조건을 벽 면에 거는 경우와, 마지막 셀 중심에 그냥 박아버리는 경우를 나란히 돌린다.

import math
 
MU, GRAD_P, RADIUS = 1.0e-3, 100.0, 1.0e-3   # Pa*s, Pa/m, m
 
 
def poiseuille_q(radius, mu=MU, grad_p=GRAD_P):
    """해석 유량  Q = pi*G*R^4/(8 mu)"""
    return math.pi * grad_p * radius ** 4 / (8.0 * mu)
 
 
def solve_pipe_fv(ncell, wall_at_face=True, radius=RADIUS, mu=MU, grad_p=GRAD_P):
    """축대칭 완전발달 속도분포, 환형 셀 ncell개의 유한체적"""
    dr = radius / ncell
    rf = [i * dr for i in range(ncell + 1)]          # 면 반지름
    lo, dg, up, rhs = [0.0] * ncell, [0.0] * ncell, [0.0] * ncell, [0.0] * ncell
    for i in range(ncell):
        rhs[i] = -grad_p * (rf[i + 1] ** 2 - rf[i] ** 2) / 2.0
        if i > 0:
            w = mu * rf[i] / dr
            lo[i], dg[i] = w, dg[i] - w
        if i < ncell - 1:
            e = mu * rf[i + 1] / dr
            up[i], dg[i] = e, dg[i] - e
        else:
            if wall_at_face:                          # 벽 면에 무활조건
                dg[i] -= mu * rf[i + 1] / (dr / 2.0)
            else:                                     # 마지막 셀 중심에 박은 무활조건
                lo[i], dg[i], up[i], rhs[i] = 0.0, 1.0, 0.0, 0.0
    for i in range(1, ncell):                         # Thomas 알고리즘
        m = lo[i] / dg[i - 1]
        dg[i] -= m * up[i - 1]
        rhs[i] -= m * rhs[i - 1]
    u = [0.0] * ncell
    u[-1] = rhs[-1] / dg[-1]
    for i in range(ncell - 2, -1, -1):
        u[i] = (rhs[i] - up[i] * u[i + 1]) / dg[i]
    q = sum(u[i] * math.pi * (rf[i + 1] ** 2 - rf[i] ** 2) for i in range(ncell))
    return q
 
 
def fit_slope(xs, ys):
    """log y 대 log x의 최소제곱 기울기"""
    lx = [math.log(x) for x in xs]
    ly = [math.log(y) for y in ys]
    n = len(lx)
    mx, my = sum(lx) / n, sum(ly) / n
    num = sum((lx[i] - mx) * (ly[i] - my) for i in range(n))
    den = sum((lx[i] - mx) ** 2 for i in range(n))
    return num / den
 
 
_seed = 20260901
 
 
def unit_normal():
    """어디서 돌려도 같은 숫자가 나오도록 LCG 위에 올린 Box-Muller"""
    global _seed
    out = []
    for _ in range(2):
        _seed = (1103515245 * _seed + 12345) % (2 ** 31)
        out.append((_seed + 0.5) / 2 ** 31)
    return math.sqrt(-2.0 * math.log(out[0])) * math.cos(2 * math.pi * out[1])
 
 
def measured_exponent(dmin, dmax, ntube, noise, repeat=200):
    """Poiseuille의 실험: 유량은 참값, 지름은 상대오차를 안고 읽는다"""
    slopes = []
    for _ in range(repeat):
        ds, qs = [], []
        for k in range(ntube):
            f = k / (ntube - 1)
            d_true = dmin * (dmax / dmin) ** f
            qs.append(poiseuille_q(d_true / 2.0))
            ds.append(d_true * (1.0 + noise * unit_normal()))
        slopes.append(fit_slope(ds, qs))
    mean = sum(slopes) / len(slopes)
    sd = math.sqrt(sum((s - mean) ** 2 for s in slopes) / len(slopes))
    return mean, sd
 
 
print("[1] wall half a cell off  (R = 1.000 mm)")
print(" N   Q_face/Q_exact   Q_centre/Q_exact   (1-1/2N)^4")
qex = poiseuille_q(RADIUS)
for n in (10, 20, 40, 80):
    a = solve_pipe_fv(n, True) / qex
    b = solve_pipe_fv(n, False) / qex
    print(f"{n:3d}      {a:.4f}            {b:.4f}          {(1-1/(2*n))**4:.4f}")
 
print()
print("[2] 1% error on D vs the resulting error on Q")
for e in (0.005, 0.01, 0.02):
    print(f"  dD/D = {e*100:4.1f}%  ->  dQ/Q = {((1+e)**4-1)*100:5.2f}%")
 
print()
print("[3] exponent fitted from 12 tubes, 200 repeats")
print(" range of D            noise on D    exponent (mean +- sd)")
for (lo_d, hi_d, tag) in ((0.10e-3, 0.30e-3, "0.10 - 0.30 mm"),
                          (0.015e-3, 0.60e-3, "0.015- 0.60 mm")):
    for nz in (0.0, 0.01, 0.03):
        m, s = measured_exponent(lo_d, hi_d, 12, nz)
        print(f" {tag}        {nz*100:4.1f}%       {m:.3f} +- {s:.3f}")

[1] wall half a cell off  (R = 1.000 mm)
 N   Q_face/Q_exact   Q_centre/Q_exact   (1-1/2N)^4
 10      1.0100            0.8100          0.8145
 20      1.0025            0.9025          0.9037
 40      1.0006            0.9506          0.9509
 80      1.0002            0.9752          0.9752
 
[2] 1% error on D vs the resulting error on Q
  dD/D =  0.5%  ->  dQ/Q =  2.02%
  dD/D =  1.0%  ->  dQ/Q =  4.06%
  dD/D =  2.0%  ->  dQ/Q =  8.24%
 
[3] exponent fitted from 12 tubes, 200 repeats
 range of D            noise on D    exponent (mean +- sd)
 0.10 - 0.30 mm         0.0%       4.000 +- 0.000
 0.10 - 0.30 mm         1.0%       3.996 +- 0.033
 0.10 - 0.30 mm         3.0%       3.993 +- 0.098
 0.015- 0.60 mm         0.0%       4.000 +- 0.000
 0.015- 0.60 mm         1.0%       4.000 +- 0.009
 0.015- 0.60 mm         3.0%       3.995 +- 0.028

첫 표의 둘째 열은 격자를 두 배로 늘릴 때마다 오차가 4분의 1로 줄어든다. 2차 수렴이다. 셋째 열은 그렇지 않다. 오차가 절반씩만 줄고, 넷째 열 $(1-1/2N)^{4}$ 과 소수점 셋째 자리까지 맞는다.

이 일치가 진단이다. 셀 중심에 무활조건을 박는 순간 계산은 반지름 $R - \Delta r/2$ 짜리 관을 푼다. 이산화 오차가 아니라 다른 관이다. 그리고 그 반 칸이 유량에서는 네 배가 된다. N = 20이면 반지름은 2.5% 작고 유량은 9.75% 작다.

벽을 노드 사이 절반에 놓는 문제는 유한체적만의 것이 아니다. 격자 볼츠만의 bounce-back도 벽을 두 노드의 중간에 세운다. 어느 쪽이든 "벽이 어디에 있다고 코드가 믿는가"를 먼저 확정해야 한다.

지수 4는 잡음 속에서도 살아남는가#

Poiseuille의 상황은 우리와 반대였다. 그는 $D$ 를 알고 $Q$ 를 예측한 것이 아니라, 잰 $Q$ 에서 지수를 읽어내야 했다. 그런데 지름은 가장 재기 어려운 양이다. 유리관 안지름 0.015 mm를 1% 정확도로 읽는 일을 생각해보면 된다.

세 번째 표가 그 실험이다. 지름에만 상대오차를 넣고 $\log Q$ – $\log D$ 기울기를 200번 다시 뽑았다. 관 지름 범위가 0.10–0.30 mm일 때 1% 잡음은 기울기에 $\pm 0.033$ 의 흔들림을 만든다. 범위를 0.015–0.60 mm로 넓히면 같은 잡음, 같은 관 개수인데 $\pm 0.009$ 로 줄어든다.

이유는 회귀식 하나로 보인다.

\sigma_{n} \approx \frac{n\,\sigma_{\delta D/D}}{\sqrt{N}\;\sigma_{\ln D}}

$n = 4$ 는 참 지수, $N$ 은 관 개수, $\sigma_{\ln D}$ 는 $\log D$ 축에서 관들이 벌어져 있는 폭이다. 좁은 범위에 넣으면 $4 \times 0.01 / (\sqrt{12} \times 0.331) = 0.035$ , 넓은 범위는 $0.010$ . 표의 0.033과 0.009가 그것이다.

지수의 정확도는 자를 정밀하게 만들어서 얻는 것이 아니라 $\log D$ 축의 지렛대를 길게 잡아서 얻는다. Poiseuille가 머리카락 5분의 1 굵기까지 내려간 이유가 여기 있다.

noise on D3.0%Dmax/Dminx3tubes12

Leave the bench running with Dmax/Dmin at 3: the histogram spreads out and the fitted exponent wanders (mean 4.000, sd 0.000 over 0 fits). Now push the range to 40 without touching the noise — the same clumsy ruler, the same number of tubes, and the histogram collapses onto 4. The exponent is bought with the lever arm in log D, not with a better ruler.

Dmax/Dmin을 3에 두고 히스토그램이 퍼지는 것을 본 다음, 잡음은 그대로 두고 범위만 40으로 밀어보자. 자도 관 개수도 그대로인데 분포가 4로 모인다.

Stokes가 $K''$ 를 $\pi/128$ 로 바꾼 뒤#

Poiseuille의 $K''$ 가 $\pi/128\mu$ 로 바뀐 것은 단위 정리가 아니었다. 실험이 남긴 빈 상수 자리에 Navier 방정식의 $\mu$ 가 들어간 사건이다. Stokes가 이 법칙을 Navier 방정식에서 수학적으로 유도해 보이면서, 관에서 잰 값과 방정식에 넣는 물성이 같은 것임이 확인됐다.

그 뒤로 이 식은 방향이 뒤집혔다. 지수를 재는 실험이 아니라, $\mu$ 를 재는 점도계가 된다. CGS 점도 단위 poise가 그의 이름에서 왔다. 방정식은 Poiseuille의 실험과 Stokes의 보강을 거쳐 Navier–Stokes 방정식으로 자리 잡았고, 2000년에는 클레이 연구소의 7대 난제 중 하나가 됐다.

우리 코드에서 같은 구조가 반복된다. 벽함수 상수, 유효 지름, 접촉각 — 실험이 형태를 주고 이론이 계수를 채운다. 형태가 $D^{4}$ 처럼 가파르면 계수보다 기하가 먼저 오차를 지배한다.

원관 검증에서 내가 먼저 보는 것#

유량이 어긋났을 때 스킴부터 열면 시간을 버린다. 순서는 이렇다.

먼저 격자를 두 배로 조밀하게 해서 오차가 4분의 1로 줄어드는지 본다. 줄지 않으면 이산화 문제가 아니다. 다음으로 오차비를 $(1-1/2N)^{4}$ 같은 기하 인자와 대조한다. 맞으면 벽 위치다. 마지막으로 유량 오차를 지름 오차로 나눠본다. 그 몫이 4에 가까우면 반지름 하나가 전부를 설명한다는 뜻이다.

Poiseuille는 지름을 알 수 없어서 지수를 읽었고, 우리는 지수를 알기 때문에 지름을 되짚을 수 있다. 같은 $D^{4}$ 를 양쪽에서 쓰는 셈이다.