Δt를 5배로 키우니 1.9배 빨라졌고, 10배에서는 제자리였다 — 암시적 표면장력이 여는 창

3일짜리 물방울 계산 앞에서#

2D 진동 물방울 하나가 사흘째 돌고 있다. 유속은 느리고 격자는 크지 않다. 그런데 시간간격이 $10^{-6}$초다. 표면장력을 명시적으로 처리하는 순간 시간간격은 유속이 아니라 모세관파가 정한다.

이럴 때 나오는 제안이 표면장력을 음함수(implicit, 새 시각의 값으로 푸는 방식)로 처리하는 것이다. 제약을 깨면 시간간격을 5배, 10배로 키울 수 있다고 한다. 그러면 사흘이 하루가 되는가.

답은 "5배까지는"이다. Janodet, van Wachem, Denner가 2025년에 낸 완전결합 알고리즘 논문은 밀도비 1000에서 이 창의 양쪽 끝을 실제로 재놓았다. 위로는 안정 한계가, 아래로는 단계당 비용이 창을 닫는다. 이 글은 그 두 벽이 각각 어디에 서 있는지, 그리고 격자를 조밀하게 해도 오차가 왜 멈추는지를 숫자로 따라간다.

시간간격을 묶는 것은 유속이 아니라 모세관파다#

계면에 표면장력이 있으면 격자가 분해할 수 있는 가장 짧은 모세관파가 하나 생긴다. 파장은 $\lambda_\sigma = 2\Delta x$다. 이 파가 한 셀을 건너가는 시간보다 큰 시간간격을 쓰면 명시적 표면장력 항은 발산한다. Denner와 van Wachem이 정리한 형태는 이렇다.

$\rho_A, \rho_B$는 두 유체의 밀도, $\sigma$는 표면장력 계수, $\Delta x$는 격자 간격이다. 지수가 $3/2$인 것이 문제다. 격자를 절반으로 줄이면 시간간격은 2.8배 작아진다. 이류 CFL 조건의 $\Delta x^1$보다 빠르게 조여든다. 확산 항은 음함수로 풀어 제약에서 빼낼 수 있지만, 표면장력은 오랫동안 그러지 못했다. 이 조건이 왜 생기고 어떻게 음함수화하는지는 모세관 시간간격 제약을 다룬 글에 따로 적었다.

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

밀도비 1000인 두 유체 사이의 모세관파가 점성으로 감쇠하는 장면이다. 회색 점선은 Prosperetti의 해석해, 파란 선은 이산 솔버가 낸 진폭이다. lambda/dx를 올리면 두 선이 붙는다. curvature dx^0.5 버튼이 켜져 있을 때 dt/dt_sigma를 0.5에서 8까지 밀어보면, 오차가 거의 움직이지 않는다는 것이 바로 다음 절의 주제다.

제약을 깬 자리에 남는 두 번째 상한#

표면장력을 음함수로 넣으면 $\Delta t_\sigma$는 넘어설 수 있다. 그런데 임의로 큰 시간간격이 되지는 않는다. Galusinski와 Vigneaux의 해석을 따라 Denner 등이 정리한 상한은 두 시간척도의 경쟁으로 나온다.

$\tau_{vc} = \hat\mu \lambda_\sigma / \sigma$는 점성-모세관 시간척도, $\tau_\sigma = \sqrt{\hat\rho \lambda_\sigma^3 / \sigma}$는 모세관 시간척도다 ($\hat\rho = \rho_A + \rho_B$, $\hat\mu = \mu_A + \mu_B$). $a_1, a_2$는 케이스마다 다른 상수로, $a_1 = 1/(16\pi)$이고 $a_2 = 0$이면 원래의 $\Delta t_\sigma$가 그대로 나온다.

두 척도의 비가 격자 온네조르게 수다.

$\mathrm{Oh}_{\Delta x} \ll 1$이면 관성이 지배해 $\Delta t^{*} \propto \tau_\sigma$, 반대면 점성이 지배해 $\Delta t^{*} \propto \tau_{vc}$다. 실무적으로 중요한 것은 후자다. 점도가 크거나 모세관파가 짧으면 상한이 크게 열린다.

숫자가 뼈아픈 쪽은 밀도비다. 정지 물방울(Laplace 평형) 케이스에서 $\mathrm{Oh}_{\Delta x} \ll 1$ 영역의 상한은 밀도비 1000일 때 $1.5\,\Delta t_\sigma$였다. 같은 알고리즘 계열이 밀도비 1에서 $15\,\Delta t_\sigma$를 썼던 것과 비교하면 열 배가 사라진 것이다. 큰 $\mathrm{Oh}_{\Delta x}$ 영역에서도 한 자릿수 차이가 난다. 현실적인 기액 밀도비는 창을 좁힌다.

격자를 8배 조밀하게 해도 오차는 절반만 줄었다#

논문의 두 번째 검증은 감쇠하는 모세관파다. 밀도비와 점도비 모두 1000, 라플라스 수 $\mathrm{La} = \rho\lambda\sigma/\mu^2 = 300$, 격자 $\lambda/\Delta x = \{25, 50, 100, 200\}$, 시간간격 $\Delta t/\Delta t_\sigma = \{0.5, 2, 8\}$. 해석해와의 차이를 진폭의 $L_2$ 노름으로 잰다.

결과 표에서 눈에 걸리는 것은 오차의 크기가 아니라 수렴차수다. 대부분의 칸이 0.46에서 0.95 사이에 있다. 같은 문제를 밀도비 1에서 풀면 2차가 나오는 자리다. 격자를 8배 조밀하게 해도 오차는 절반쯤만 줄어든다.

논문이 지목한 원인은 시간 이산화가 아니다. 계면 수송이다. 논리는 두 줄이다. 쓰이는 계면 포착 스킴은 잘해야 2차 정확도다. 곡률은 색함수의 2계 미분이므로 정확도가 두 단계 떨어진다. 즉 곡률은 잘해야 0차다. 충분히 조밀한 격자에서 진폭 오차의 수렴차수는 결국 0으로 간다. 오차가 상수로 눌러앉고 더 이상 줄지 않는다는 뜻이다.

Python으로 확인한 수렴차수의 주인#

이 논리를 감쇠 진동자 하나로 줄여 확인할 수 있다. 모세관파의 진폭은 선형 영역에서 $A'' + 2\nu k^2 A' + \omega_0^2 A = 0$을 따른다. 솔버가 보는 것은 $\omega_0$가 아니라 곡률 오차가 섞인 $\omega_{num} = \omega_0\sqrt{1 + C(\Delta x/\lambda)^q}$다. $q$만 바꿔가며 사다리꼴 적분으로 마치고, 해석해와의 $L_2$ 노름과 수렴차수를 뽑는다.

import math
 
SIGMA, RHO_HAT, LAMBDA, K, LA = 1.0, 1.0, 2*math.pi, 1.0, 300.0
MU = math.sqrt(RHO_HAT * LAMBDA * SIGMA / LA)
NU = MU / RHO_HAT
A0, T_END = LAMBDA / 100.0, 25.0
 
 
def capillary_omega(dx, q, c_kappa=0.6):
    """이산 솔버가 실제로 보는 진동수 — 곡률 오차는 O(dx^q)"""
    w0 = math.sqrt(SIGMA * K**3 / RHO_HAT)
    return w0 * math.sqrt(1.0 + c_kappa * (dx / LAMBDA) ** q)
 
 
def analytic_amplitude(t):
    """A'' + 2*nu*k^2*A' + w0^2*A = 0 의 해석해"""
    w0 = math.sqrt(SIGMA * K**3 / RHO_HAT)
    g = NU * K**2
    wd = math.sqrt(w0**2 - g**2)
    return A0 * math.exp(-g*t) * (math.cos(wd*t) + g/wd * math.sin(wd*t))
 
 
def march_amplitude(dt, w, n_steps):
    """[A, A'] 를 사다리꼴(Crank-Nicolson)로 전진"""
    g = NU * K**2
    a, v, hist = A0, 0.0, [A0]
    for _ in range(n_steps):
        h = 0.5 * dt
        rhs_a, rhs_v = a + h*v, v + h*(-w**2 * a - 2*g*v)
        det = (1 + 2*g*h) + h*h*w**2
        a = ((1 + 2*g*h) * rhs_a + h * rhs_v) / det
        v = (-h * w**2 * rhs_a + rhs_v) / det
        hist.append(a)
    return hist
 
 
def l2_amplitude(hist, dt):
    """진폭의 L2 오차 노름 (논문 Eq. 61)"""
    acc = 0.0
    for i, a in enumerate(hist):
        w = 0.5 if i in (0, len(hist)-1) else 1.0
        acc += w * (a - analytic_amplitude(i*dt))**2 * dt
    return math.sqrt(acc / (len(hist)-1) / dt) / A0
 
 
def order_of(e_coarse, e_fine):
    return math.log(e_coarse / e_fine) / math.log(2.0)
 
 
for label, q in [("curvature error ~ dx^2", 2.0), ("curvature error ~ dx^0.5", 0.5)]:
    print(f"\n{label}")
    print("lam/dx |  dt/dt_s=0.5        dt/dt_s=2          dt/dt_s=8")
    prev = {}
    for n in [25, 50, 100, 200]:
        dx = LAMBDA / n
        dt_sigma = math.sqrt(RHO_HAT * dx**3 / (2*math.pi*SIGMA))
        w = capillary_omega(dx, q)
        row = []
        for s in [0.5, 2.0, 8.0]:
            dt = s * dt_sigma
            e = l2_amplitude(march_amplitude(dt, w, int(T_END/dt)), dt)
            tag = "  (-- )" if s not in prev else f" ({order_of(prev[s], e):4.2f})"
            row.append(f"{e:.3e}{tag}")
            prev[s] = e
        print(f"{n:6d} | " + "  ".join(row))

curvature error ~ dx^2
lam/dx |  dt/dt_s=0.5        dt/dt_s=2          dt/dt_s=8
    25 | 5.763e-04  (-- )  5.692e-04  (-- )  1.680e-02  (-- )
    50 | 1.516e-04 (1.93)  6.731e-05 (3.08)  2.023e-03 (3.05)
   100 | 3.885e-05 (1.96)  2.548e-05 (1.40)  2.345e-04 (3.11)
   200 | 9.833e-06 (1.98)  8.085e-06 (1.66)  2.506e-05 (3.23)
 
curvature error ~ dx^0.5
lam/dx |  dt/dt_s=0.5        dt/dt_s=2          dt/dt_s=8
    25 | 7.565e-02  (-- )  7.483e-02  (-- )  6.051e-02  (-- )
    50 | 5.449e-02 (0.47)  5.439e-02 (0.46)  5.266e-02 (0.20)
   100 | 3.897e-02 (0.48)  3.896e-02 (0.48)  3.874e-02 (0.44)
   200 | 2.775e-02 (0.49)  2.775e-02 (0.49)  2.773e-02 (0.48)

윗 표에서는 시간간격이 작을 때 2차가 나온다. $\Delta t/\Delta t_\sigma = 8$ 열의 3차는 $\Delta t \propto \Delta x^{3/2}$이므로 2차 시간오차가 $\Delta x^3$으로 줄기 때문이다.

아랫 표가 논문이 말한 상황이다. 수렴차수가 0.5 부근에 고정된다. 더 중요한 것은 세 열의 값이 사실상 같다는 점이다. 시간간격을 16배 줄여도 오차가 그대로다. 정확도의 바닥을 정하는 것이 시간 이산화가 아니라 곡률이라는 뜻이다. 논문의 실측 수렴차수 0.46~0.95는 이 그림과 정확히 겹친다.

이 지점은 계면 이류의 CFL 상한을 다룬 글과 짝을 이룬다. 거기서는 시간간격을 키워도 CFL 0.05가 남았고, 여기서는 격자를 줄여도 곡률이 남는다.

시간간격을 5배로 키우면 1.9배, 10배로 키우면 그대로#

세 번째 검증은 2D 타원 물방울의 감쇠 진동이다. 장반경 0.15 m, 단반경 0.1 m에서 출발해 $n=2$ 모드로 진동하다 점성으로 멎는다. 적용 시간간격은 두 제약의 최솟값이다.

$\Sigma$는 모세관 제약을 몇 배로 깰지 정하는 인자다. 논문은 $\Sigma \in \{2, 5, 10\}$을 돌렸고 CFL 상한은 0.05로 묶었다.

정확도부터 보면, 진동 주파수 오차는 $\Sigma = 2$와 $5$에서 약 3%였다. 같은 해상도에서 명시적 표면장력으로 얻은 약 4.5%보다 작다. 반면 $\Sigma = 10$은 운동에너지 감쇠를 제대로 따라가지 못했다. 표면장력이 이끄는 계면 운동이 시간적으로 분해되지 않으니, 시간 이산화의 형식적 2차 정확도를 기대할 수 없다는 것이 논문의 설명이다.

비용 쪽 숫자가 이 글의 제목이다. $\Sigma$를 2에서 5로, 즉 2.5배 올렸을 때 전체 벽시계 시간은 1.9배 줄었다. 그런데 10으로 더 올렸을 때 이 이득은 유지되지 않았다. 단계당 환산 계산시간이 크게 늘었기 때문이다. 원인은 하나다. 시간간격이 커질수록 매 스텝의 비선형 반복이 느리게 수렴한다. 단계 수는 줄지만 단계당 비용이 그만큼 오른다.

S 슬라이더가 $\Sigma$다. 1에서 5까지 올리면 파란 레인이 눈에 띄게 먼저 끝나지만, 10까지 밀면 거의 제자리다. Oh_dx를 내리면 붉은 벽($\Delta t^{*}$)이 왼쪽으로 와서 빠른 레인이 아예 죽는다.

벽들은 서로 다른 곳에 서 있다#

이 케이스 하나에 걸린 제약과 최적점을 한 표에 세우면 다섯 줄이다. 성격이 전부 다르다.

$\Delta t_\sigma$를 깨는 알고리즘은 이 표의 두 번째 줄만 지운다. 나머지는 그대로 남는다. 정지 물방울에서 힘 균형이 어떻게 유지되는지는 기생 전류를 다룬 글에서 다뤘다.

그래서 $\Sigma$는 어떻게 고르는가#

논문의 결론은 $\Sigma$에 케이스마다 다른 최적값이 있다는 것이고, 이 문제에서는 5였다. 그 값을 찾는 순서는 셋이다.

먼저 $\mathrm{Oh}_{\Delta x}$를 계산한다. 이 수가 1보다 훨씬 작으면 안정 창 자체가 좁다. 밀도비가 크면 더 좁다. $\Sigma$를 10으로 놓고 시작할 이유가 없다.

다음으로 물리 시간척도를 센다. 관심 있는 진동 모드의 주기를 몇 스텝으로 자르고 있는지 세어보면 된다. 안정하다고 정확한 것이 아니다. $\Sigma = 10$은 안정했지만 에너지 감쇠를 놓쳤다.

마지막으로 스텝당 비선형 반복 횟수를 로그에서 읽는다. $\Sigma$를 올렸는데 반복 횟수가 비례해서 늘면 그 지점이 창의 오른쪽 끝이다. 벽시계 시간은 이미 바닥을 지났다.