점성이 소용돌이를 죽이지 않는다 — 뉴턴의 1/r과 Lamb–Oseen 코어

점성 솔버를 검증할 때 흔히 쓰는 초기조건이 자유 와류다. 속도장 $v_\theta = \Gamma/2\pi r$를 넣고 시간을 돌린다. 점성이 있으니 와류는 천천히 약해질 것이다. 그런데 반지름 0.5 지점의 속도가 50배의 시간이 지나도 소수점 넷째 자리까지 그대로다. 코드가 틀린 게 아니다. 이 글은 왜 그런지, 그리고 실제로 와류를 죽이는 것이 무엇인지를 다룬다. 답은 1687년에 이미 절반이 나와 있었다.

1687년, 프린키피아 2권이 겨눈 것#

『프린키피아』는 세 권이다. 1권은 힘의 법칙, 3권은 만유인력이다. 물리학자들이 거의 언급하지 않는 2권이 유체역학이다. 표적이 분명했다. 데카르트의 소용돌이(vortex) 우주론이다.

데카르트는 행성이 태양 주위를 도는 것을 유체 소용돌이에 실려 도는 것으로 설명했다. 우주는 미세한 물질로 가득 차 있고, 그 물질이 거대한 회전 흐름을 이룬다는 그림이다. 뉴턴은 이 그림을 유체역학으로 반박했다. 유체의 운동에는 반드시 저항(resistance)이 따르므로, 외력이 없으면 소용돌이는 결국 소멸한다는 논지였다.

논지를 세우려면 저항을 정량화해야 했다. 그래서 2권에서 뉴턴은 유체의 저항이 전단율에 선형으로 비례한다고 놓는다.

$\tau$는 전단응력, $du/dy$는 속도 구배, 비례계수 $\mu$가 점도다. 유체 점도가 수학적으로 정의된 최초의 자리다. 후대에 이 선형 관계를 따르지 않는 유체가 훨씬 많다는 게 밝혀지면서, 따르는 쪽을 뉴턴 유체라 부르게 됐다.

실린더 하나에서 나오는 1/r#

뉴턴의 유도는 지금 봐도 현대적이다. 무한히 긴 실린더가 점성 유체 안에서 일정한 각속도로 돈다. 정상 상태에서 반지름 $r$인 원통면을 통과해 바깥으로 전달되는 토크는 어디서나 같아야 한다. 안 그러면 그 사이 층에 각운동량이 쌓인다.

원통 좌표계에서 순수 회전 유동의 전단응력은 $\tau_{r\theta} = \mu\, r\, \frac{d}{dr}(v_\theta/r)$다. 단순히 $\mu\, dv_\theta/dr$가 아니다. 강체 회전($v_\theta \propto r$)에는 전단이 없어야 하는데, 이 항이 그 조건을 자동으로 만족시킨다. 단위 길이당 토크는 응력에 팔 $r$과 둘레 $2\pi r$을 곱한 값이다.

이 상미분방정식을 풀면 두 항이 나온다.

첫 항은 강체 회전, 둘째 항은 자유 와류다. 바깥 경계가 무한히 멀고 그곳에서 유체가 정지해 있으면 $A=0$이다. 남는 것은 $v_\theta \propto 1/r$이다. 이것이 뉴턴이 2권에서 얻은 결과이고, 오늘날 Taylor–Couette 유동의 정상해로 배우는 바로 그 식이다.

케플러가 요구하는 지수는 1/2다#

이제 반박이 성립한다. 데카르트의 소용돌이가 행성을 나른다면, 그 소용돌이의 속도 분포는 케플러 제3법칙을 재현해야 한다. 공전 주기가 $T \propto r^{3/2}$이므로 속도는 $v = 2\pi r / T \propto r^{-1/2}$다.

뉴턴의 유체 소용돌이는 $r^{-1}$을 준다. 케플러는 $r^{-1/2}$을 요구한다. 두 지수가 다르다. 점성 유체의 정상 소용돌이는 행성 궤도를 만들 수 없다.

아래 다이얼에서 지수를 직접 움직여 보자.

n = −1에서 스포크가 곧게 유지되는 것이 강체 회전이다. n = 1로 밀면 스포크가 감기는데도 오른쪽 아래 "viscous force" 막대는 초록으로 떨어진다. 케플러의 n = 0.5에서는 두 막대가 모두 빨갛다. 관찰 포인트는 이 지점이다 — 그 프로파일을 유지하려면 누군가 계속 힘을 넣어 줘야 한다.

응력은 있는데 힘이 없다#

여기서 처음의 의문이 풀린다. $v_\theta = C\,r^{-n}$을 회전 성분의 점성 항에 넣어 보자.

마지막 $-v_\theta/r^2$ 항은 곡률 때문에 붙는다. 직교좌표의 라플라시안에는 없는 항이다. 괄호를 정리하면 $n^2-1$이 남는다.

이 계수는 $n = \pm 1$에서 정확히 0이다. $n = -1$은 강체 회전이다. $n = +1$이 자유 와류다. 즉 자유 와류에는 점성력이 한 점도 작용하지 않는다.

전단응력이 0이라서가 아니다. 응력 자체는 $\tau_{r\theta} = -\mu\Gamma/\pi r^{2}$로 코어 근처에서 대단히 크다. 다만 한 유체 요소의 안쪽 면이 받는 토크와 바깥쪽 면이 받는 토크가 정확히 상쇄된다. 힘은 응력이 아니라 응력의 발산이다. 자유 와류는 그 발산이 0인 특수한 프로파일이다.

같은 이야기를 와도로도 할 수 있다. $\omega = \frac{1}{r}\frac{d(r v_\theta)}{dr}$이고, $v_\theta = \Gamma/2\pi r$이면 $r v_\theta$가 상수라서 $r>0$ 어디서나 $\omega = 0$이다. 비압축성 유동의 점성력은 $-\nu\,\nabla\times\boldsymbol{\omega}$로 쓸 수 있다. 와도가 없으면 점성력도 없다.

코어가 퍼지는 속도로 시간이 정해진다#

와도가 없다고 했지만 정확히는 원점을 제외한 곳이다. 순환 $\Gamma$는 어딘가에 있어야 하고, 이상적인 자유 와류에서는 그것이 원점의 델타 함수로 몰려 있다. 점성이 실제로 일하는 자리가 여기다.

원점에 순환 $\Gamma$를 집중시켜 놓고 확산 방정식을 풀면 Lamb–Oseen 와류가 나온다.

지수 항이 코어를 만든다. $r \gg \sqrt{4\nu t}$에서는 괄호가 1이 되어 자유 와류로 돌아간다. 코어 반경은 $r_c = 2.2418\sqrt{\nu t}$로 자라고, 최대 속도는 $t^{-1/2}$로 떨어진다.

핵심은 순환이다. $r \to \infty$에서 순환은 언제나 정확히 $\Gamma$다. 점성은 와도를 퍼뜨릴 뿐 없애지 않는다. 아래 실험에서 벽 스위치를 켜고 꺼 보자.

벽이 꺼져 있으면 코어는 부풀고 피크는 내려가지만 오른쪽 아래 순환 막대는 1.000에서 움직이지 않는다. 벽을 켜면 같은 막대가 내려가기 시작한다. 조작할 것은 $\nu$와 벽, 관찰할 것은 코어 반경과 순환 막대의 차이다.

Python으로 세어본 반경별 감쇠율#

말로 한 것을 숫자로 확인하자. 축대칭 회전 방정식 $\partial_t v_\theta = \nu(\partial_{rr} v_\theta + \partial_r v_\theta / r - v_\theta / r^2)$를 셀 중심 격자 400개에서 전진 오일러로 적분한다.

import numpy as np
 
NU, GAMMA, R, N = 1.0e-3, 1.0, 1.0, 400
dr = R / N
r = (np.arange(N) + 0.5) * dr
 
def lamb_oseen(rr, t):
    """Lamb-Oseen 와류의 접선 속도. t가 작으면 맨 1/r 와류에 수렴한다."""
    return GAMMA / (2 * np.pi * rr) * (1 - np.exp(-rr * rr / (4 * NU * t)))
 
def swirl_terms(v):
    """회전 라플라시안의 세 조각. 상쇄를 눈으로 보려고 따로 둔다."""
    ghost = np.concatenate(([-v[0]], v, [v[-1] * r[-1] / (r[-1] + dr)]))
    return ((ghost[2:] - 2 * ghost[1:-1] + ghost[:-2]) / dr**2,
            (ghost[2:] - ghost[:-2]) / (2 * dr) / r,
            -v / r**2)
 
def swirl_operator(v, outer):
    """nu * (v_rr + v_r/r - v/r^2): 순수 회전에 남는 점성 항 전부."""
    g_out = -v[-1] if outer == 'wall' else v[-1] * r[-1] / (r[-1] + dr)
    w = np.concatenate(([-v[0]], v, [g_out]))
    v_rr = (w[2:] - 2 * w[1:-1] + w[:-2]) / dr**2
    v_r = (w[2:] - w[:-2]) / (2 * dr)
    return NU * (v_rr + v_r / r - v / r**2)
 
def march_swirl(v, t_end, outer):
    dt = 0.2 * dr * dr / NU
    for _ in range(int(round(t_end / dt))):
        v = v + dt * swirl_operator(v, outer)
    return v
 
def circulation_at(v, radius):
    return 2 * np.pi * radius * np.interp(radius, r, v)
 
# 1. 맨 1/r 와류: 큰 항 세 개가 서로 지운다
a, b, c = (np.interp(0.10, r, x) for x in swirl_terms(GAMMA / (2 * np.pi * r)))
print(f"free vortex at r=0.10:  v_rr={a:+8.2f}  v_r/r={b:+8.2f}  -v/r^2={c:+8.2f}  sum={a+b+c:+.2e}")
print(f"                        shear stress tau_rtheta = {-NU * GAMMA / (np.pi * 0.10**2):+.4f}")
 
# 2. 실제로 전진시켜서 해석해와 대조
t0, t1 = 0.02, 1.0
v0 = lamb_oseen(r, t0)
v1 = march_swirl(v0, t1 - t0, 'free')
print(f"\nmarched {t0} -> {t1} s   L-inf vs Lamb-Oseen = {np.max(np.abs(v1 - lamb_oseen(r, t1))):.1e}")
print("    r    v(0.02)   v(1.00)    change")
for x in (0.02, 0.05, 0.15, 0.50, 0.95):
    p, q = np.interp(x, r, v0), np.interp(x, r, v1)
    print(f"{x:5.2f} {p:9.4f} {q:9.4f} {100 * (q - p) / p:+8.1f}%")
 
# 3. 코어는 sqrt(nu t)로 퍼지고 피크는 1/sqrt(t)로 떨어진다. 순환은 그대로다
print("\n    t   r_core   v_peak  v_peak*sqrt(t)  Gamma(0.9)")
v, tc = v0.copy(), t0
for t in (0.05, 0.20, 0.50, 1.00):
    v = march_swirl(v, t - tc, 'free'); tc = t
    i = int(np.argmax(v))
    print(f"{t:5.2f} {r[i]:8.4f} {v[i]:8.4f} {v[i] * np.sqrt(t):13.4f} {circulation_at(v, 0.9):11.4f}")
 
# 4. r = R 에 벽을 세우면 같은 와류가 죽는다
print("\n    t   Gamma(0.9)  unbounded    walled")
vw, tc = v0.copy(), t0
for t in (1.0, 20.0, 70.0, 200.0):
    vw = march_swirl(vw, t - tc, 'wall'); tc = t
    print(f"{t:8.1f}     {circulation_at(lamb_oseen(r, t), 0.9):9.4f} {circulation_at(vw, 0.9):9.4f}")
print(f"\nslowest walled mode:  R^2/(nu*j11^2) = {R**2 / (NU * 3.8317**2):.1f} s")

free vortex at r=0.10:  v_rr= +318.81  v_r/r= -159.40  -v/r^2= -159.30  sum=+9.98e-02
                        shear stress tau_rtheta = -0.0318
 
marched 0.02 -> 1.0 s   L-inf vs Lamb-Oseen = 7.7e-04
    r    v(0.02)   v(1.00)    change
 0.02    7.9233    0.7562    -90.5%
 0.05    3.1851    1.4782    -53.6%
 0.15    1.0611    1.0573     -0.4%
 0.50    0.3183    0.3183     +0.0%
 0.95    0.1675    0.1675     +0.0%
 
    t   r_core   v_peak  v_peak*sqrt(t)  Gamma(0.9)
 0.05   0.0163   7.1699        1.6032      1.0000
 0.20   0.0312   3.5880        1.6046      1.0000
 0.50   0.0513   2.2697        1.6049      1.0000
 1.00   0.0713   1.6057        1.6057      1.0000
 
    t   Gamma(0.9)  unbounded    walled
     1.0        1.0000    0.9773
    20.0        1.0000    0.4182
    70.0        0.9446    0.2186
   200.0        0.6367    0.0342
 
slowest walled mode:  R^2/(nu*j11^2) = 68.1 s

첫 블록이 상쇄를 보여준다. 세 항이 각각 300 규모인데 합은 0.1이다. 상대적으로 $3\times10^{-4}$이고, 이건 2차 차분의 이산화 오차다. 해석적으로는 정확히 0이다. 같은 자리의 전단응력은 0이 아니다.

둘째 블록이 처음의 의문에 답한다. $r = 0.5$와 $r = 0.95$에서 변화가 0.0%다. 50배의 시간을 돌려도 그렇다. $r = 0.02$에서만 90% 깎였다. 점성이 일한 곳은 코어뿐이다.

셋째 블록의 마지막 두 열이 스케일링을 확인해 준다. v_peak*sqrt(t)가 1.603에서 1.606까지 0.2% 안에서 고정된다. 코어 반경도 $2.2418\sqrt{\nu t} = 0.0709$와 0.0713으로 맞는다. 그동안 순환은 1.0000이다.

벽을 세우면 같은 와류가 죽는다#

넷째 블록이 결론이다. 무한 영역에서는 $t = 200$에서도 $r=0.9$의 순환이 0.6367로 남아 있다. 이건 감쇠가 아니라 코어가 그 반경까지 부푼 결과다. 반경을 더 키워 재면 여전히 1이다.

벽을 세운 쪽은 0.0342다. 95% 이상이 사라졌다. 벽이 하는 일은 각운동량을 계 밖으로 빼내는 것이다. 감쇠 시간은 가장 느린 모드가 정한다. $v_\theta = 0$과 축 정칙성을 만족하는 모드는 $J_1(j_{1,1} r/R)$이고, 시간 상수는 $R^2/(\nu\, j_{1,1}^2) = 68.1$초다. 계산에서 $t = 70$일 때 0.2186까지 떨어졌으니 대략 맞는다.

그러니 뉴턴의 문장은 이렇게 다듬어야 한다. 소용돌이를 죽이는 것은 점성 자체가 아니다. 점성이 있고 각운동량이 빠져나갈 경계가 있을 때 소용돌이가 죽는다. 무한 영역에서 점성은 와도를 퍼뜨리기만 한다. 데카르트를 반박하는 데는 그걸로 충분했다. 우주가 유한하든 무한하든, 소용돌이 안의 행성은 케플러의 지수를 만들 수 없다.

와류 테스트 케이스를 짜기 전에#

이 결과는 실무에서 세 가지로 쓰인다.

정확도 검증에 자유 와류를 쓰지 말 것. 코어 바깥은 점성 항이 정확히 0이라, 그 영역의 오차는 대류 스킴의 오차이지 점성 이산화의 오차가 아니다. 점성 항을 검증하려면 Lamb–Oseen의 코어를 격자로 충분히 덮고 $r_c(t)$의 성장을 재야 한다.

수치 감쇠를 여기서 잰다. 반대로 이 성질은 좋은 진단이 된다. 자유 와류를 넣고 코어 밖에서 속도가 떨어지면, 그건 물리가 아니라 코드의 수치 소산이다. 상류 차분 계열은 여기서 바로 들킨다.

계산 영역의 크기가 답을 바꾼다. 와류가 오래 남아야 하는 문제 — 날개 끝 와류, 로터 후류 — 에서 바깥 경계를 가까이 두면 벽처럼 작동한다. $R^2/(\nu\, j_{1,1}^2)$을 물리 시간과 비교해 보고 영역을 잡는 편이 낫다. 난류 점성을 쓰는 계산이면 $\nu$ 자리에 $\nu_t$를 넣어야 하고, 그러면 이 시간 상수가 수백 배 짧아진다.