Δ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$보다 빠르게 조여든다. 확산 항은 음함수로 풀어 제약에서 빼낼 수 있지만, 표면장력은 오랫동안 그러지 못했다. 이 조건이 왜 생기고 어떻게 음함수화하는지는 모세관 시간간격 제약을 다룬 글에 따로 적었다.

dt/dt_sigma를 0.5에서 8로 밀면, 곡률 오차가 켜져 있을 때 파란 진폭이 회색 점선에서 거의 안 벗어난다.

밀도비 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}$ 영역에서도 한 자릿수 차이가 난다. 현실적인 기액 밀도비는 창을 좁힌다.

격자를 한 번 반으로 줄이면 $\Delta t_\sigma$는 $2^{3/2}\approx 2.83$배가 된다. $\lambda/\Delta x$를 25에서 200으로 올리는 일은 간격을 8배 줄이는 일이다. $8^{3/2}\approx 22.6$이므로, 같은 물리 시간을 명시적 표면장력으로 밀려면 스텝 수가 약 23배가 된다. 사흘이 두 달을 넘긴다.

음함수로 바꿔도 밀도비 1000의 정지 물방울은 상한이 $1.5\,\Delta t_\sigma$에서 닫힌다. 밀도비 1에서 같은 계열이 쓰던 $15\,\Delta t_\sigma$의 10분의 1이다. $\Sigma=5$는 그 상한 밖이다. 관성 지배에서는 5배도 발산이고, 10배는 더 바깥이다. 점도가 커서 $\Delta t^{*}\propto\tau_{vc}$로 넘어가야 5배가 안정 영역에 들어온다. $a_2=0$이고 $a_1=1/(16\pi)$이면 $\Delta t^{*}$는 원래의 $\Delta t_\sigma$로 돌아간다. 창을 여는 항은 $a_2\tau_{vc}$다.

격자를 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)

S를 1에서 5로 올리면 파란 레인이 먼저 끝나고, 10으로 밀면 거의 안 움직인다.

Oh_dx를 내리면 붉은 벽 $\Delta t^{*}$가 왼쪽으로 온다. 빠른 레인이 안정 상한에 걸려 멈춘다. 아랫표는 시간간격을 16배 벌려도 오차가 그대로인 다른 벽이다.

윗 표에서는 시간간격이 작을 때 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는 이 그림과 정확히 겹친다.

lam/dx = 200 행은 2.775e-02, 2.775e-02, 2.773e-02다. dt/dt_s를 0.5에서 8로 열여섯 배 벌려도 끝자리만 움직이고, 차수는 0.49, 0.49, 0.48이다. 0.46에서 0.95 사이에 그대로 들어간다. 윗표 dt/dt_s = 0.5 열은 5.763e-04에서 9.833e-06까지 떨어지고 차수가 1.93, 1.96, 1.98로 2차에 붙는다. 곡률 오차의 지수 q를 2에서 0.5로 내리는 것만으로 바닥이 바뀐다. 시간간격 열을 갈아 생기는 차이가 아니다.

그 곡률을 체적힘으로 넣는 정식은 접촉각과 CSF 쪽에서 이어진다. 격자를 조여도 곡률의 0차가 남는다.

시간간격을 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으로 더 올렸을 때 이 이득은 유지되지 않았다. 단계당 환산 계산시간이 크게 늘었기 때문이다. 원인은 하나다. 시간간격이 커질수록 매 스텝의 비선형 반복이 느리게 수렴한다. 단계 수는 줄지만 단계당 비용이 그만큼 오른다. 2.5배를 올려 벽시계가 1.9배만 준 나머지 0.6배는, 스텝마다 늘어난 뉴턴 반복이 먹은 것이다.

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

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

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

$\Sigma=10$은 안정해도 감쇠를 놓친다#

$\mathrm{Oh}_{\Delta x}$가 1보다 훨씬 작으면 안정 창은 $1.5\,\Delta t_\sigma$ 근처에서 닫힌다. 밀도비 1000의 정지 물방울이 그 경우다. 여기서 $\Sigma$를 10으로 시작하면 창 밖이다. 진동 한 주기를 몇 스텝으로 자르는지도 같이 봐야 한다. $\Sigma=2$와 $5$의 주파수 오차는 약 3%로, 명시적 처리의 약 4.5%보다 작았다. $\Sigma=10$은 발산하지 않았는데도 운동에너지 감쇠를 따라가지 못했다.

로그의 비선형 반복이 $\Sigma$에 비례해 늘면 벽시계는 이미 바닥을 지난 것이다. 2에서 5로, 2.5배를 올릴 때 시간은 1.9배 줄었고, 10에서는 그 이득이 되돌아갔다. 단계 수는 줄었지만 단계당 뉴턴 반복이 빈자리를 채웠다. CFL은 이 계산에서 0.05로 고정돼 있다. $\Sigma$를 더 올려도 이류 쪽이 먼저 시간간격을 막는 격자가 있다.