물방울 시뮬레이션이 Δt에 발목 잡히는 이유 — 표면장력을 음함수로 푼 논문

1.5배 스텝에서 명시적 피크는 3.6e+07, 음함수는 1.000이다. 족쇄를 푸는 쪽은 후방 오일러다.

격자를 두 배로 촘촘히 하면 시간 간격은 약 2.8배가 줄어든다. 물리가 빨라진 게 아니다. 표면장력이 있는 계면에서는 모세관 시간 간격(capillary time-step constraint)이 거의 항상 제일 빡센 제약이다.

$\Delta t/\Delta t_\sigma$ 슬라이더를 0.9 위로 밀면 주황(명시적)만 진폭이 터진다. 시안(음함수)은 같은 정상파를 붙들고 있다.

Δt / Δt_σ = 0.60mesh Oh = 0.10

explicit surface tension: stableimplicit: stable—— analytic decay

점선이 해석해다. mesh Oh를 올리면 감쇠가 세져 둘 다 더 빨리 잦아든다. Janodet·van Wachem·Denner의 2025년 결합은 이 주황이 터지는 지점을 행렬 안으로 넣는다.

리뷰 대상은 R. Janodet, B. van Wachem, F. Denner, "A fully-coupled algorithm with implicit surface tension treatment for interfacial flows with large density ratios", arXiv:2410.17757 (2024).

모세관 시간 간격이라는 족쇄#

비압축성 계면 유동 솔버는 보통 확산 항을 음함수로 풀어 점성 시간 제약을 지운다. 그러면 남는 가장 강한 제약이 표면장력이다. Brackbill이 1992년에 처음 쓴 형태로,

\Delta t_\sigma = \sqrt{\frac{(\rho_A + \rho_B)\,\Delta x^3}{2\pi\sigma}} \propto \Delta x^{3/2}

$\rho_A, \rho_B$ 는 두 유체의 밀도, $\sigma$ 는 표면장력 계수, $\Delta x$ 는 격자 간격이다. 지수 $3/2$ 가 문제다. 계면을 또렷하게 보려고 격자를 절반으로 줄이면 시간 간격은 $2^{1.5}\approx 2.8$ 배 줄어든다. 작은 물방울, 모세관 규모 유동일수록 스텝 수가 폭발한다. 명시적 표면장력으로는 빠져나갈 길이 없다.

왜 Δx^1.5인가 — 모세관파 분산관계#

이 제약은 사실 CFL 조건이다. 계면 위 작은 정현파 교란은 모세관파(capillary wave, 표면장력이 복원력인 파동)로 진동한다. 비점성 심해 극한에서 분산관계는

\omega^2 = \frac{\sigma k^3}{\rho_A + \rho_B}

$\omega$ 는 각진동수, $k$ 는 파수다. 격자가 분해할 수 있는 가장 짧은 파는 $k_{\max} = \pi/\Delta x$ 다. 이 파의 진동 주기보다 $\Delta t$ 가 크면 명시적 적분이 그 진동을 따라잡지 못하고 발산한다. $\omega_{\max}\Delta t_\sigma$ 를 계산하면 $\sqrt{(k_{\max}\Delta x)^3/(2\pi)} = \pi/\sqrt{2}\approx 2.22$ 로 떨어진다. 즉 $\Delta t_\sigma$ 는 가장 빠른 모세관파에 걸린 CFL 그 자체다.

표면장력을 음함수로 — 한 줄짜리 진동자로 본 핵심#

복잡한 솔버를 다 들어낼 필요는 없다. 모세관파 한 모드는 감쇠 조화진동자로 환원된다.

\ddot{\eta} + 2\gamma\dot{\eta} + \omega^2 \eta = 0

$\eta$ 는 계면 진폭, $\gamma$ 는 점성 감쇠다. 명시적 표면장력은 이 진동자를 심플렉틱 오일러로 푸는 것과 같고, $\omega\Delta t > 2$ 에서 발산한다 — 바로 위의 모세관 CFL이다. 표면장력 항을 새 시각의 값으로 두면(후방 오일러) 무조건 안정해진다.

앞의 주황이 0.9 근처에서 터지는 이유는 심플렉틱 오일러의 한계 $\omega\Delta t = 2$ 를, 출력에 찍히는 $\omega\Delta t_\sigma = 2.221$ 로 나눈 값이기 때문이다. $2/2.221 \approx 0.900$ . 그 위로는 명시적이 진폭을 키운다. 시안은 3배까지 점선을 따라 감쇠한다. mesh Oh(점성 척도)를 키우면 감쇠가 세져 둘 다 빨리 잦아든다.

정지한 계면에 남는 가짜 흐름은 이 복원력의 이산화가 어긋날 때 생긴다. 힘 균형이 깨진 기생 전류가 그 자리다.

증폭인자로 본 안정성 경계#

진동자의 증폭인자(amplification factor, 한 스텝 뒤 진폭이 곱해지는 배율) $|G|$ 로 경계가 또렷해진다. numpy 30줄로 재현한다. 토이 문제는 충격파 튜브가 아니라 정상 모세관파의 감쇠다.

import numpy as np
 
def capillary_omega(sigma, k, rho_hat):
    # 분산관계: 비점성 심해 모세관파 omega^2 = sigma k^3 / (rhoA+rhoB)
    return np.sqrt(sigma * k**3 / rho_hat)
 
def dt_capillary(sigma, dx, rho_hat):
    # 식(4): Brackbill 모세관 시간 간격 ~ dx^{3/2}
    return np.sqrt(rho_hat * dx**3 / (2*np.pi*sigma))
 
def run_standing_wave(omega, gamma, dt, steps, implicit):
    eta, v = 1.0, 0.0
    peak = 0.0
    for _ in range(steps):
        if implicit:        # 후방 오일러: 표면장력 항을 새 시각 값으로 (무조건 안정)
            v = (v - dt*omega**2*eta) / (1 + dt**2*omega**2 + 2*gamma*dt)
            eta = eta + dt*v
        else:               # 심플렉틱 오일러: 명시적 표면장력 (omega*dt<2 제약)
            v = v - dt*omega**2*eta - dt*2*gamma*v
            eta = eta + dt*v
        peak = max(peak, abs(eta))
    return peak
 
sigma, rho_hat, L = 0.07, 1001.2, 1.0      # 물/공기
for N in (64, 128, 256):
    dx = L/N; k = np.pi/dx                  # 가장 짧은 분해 가능한 모세관파
    om = capillary_omega(sigma, k, rho_hat)
    dts = dt_capillary(sigma, dx, rho_hat)
    print(f"N={N:4d}  dt_sigma={dts:.2e}  omega*dt_sigma={om*dts:.3f}")
# 출력: 모든 N에서 omega*dt_sigma = 2.221 (= pi/sqrt(2)), CFL이 곧 모세관 제약
 
N = 128; dx = L/N; k = np.pi/dx
om = capillary_omega(sigma, k, rho_hat); dts = dt_capillary(sigma, dx, rho_hat)
dt = 1.5*dts; gamma = 0.1*om                # 모세관 제약의 1.5배로 강제
pe = run_standing_wave(om, gamma, dt, 200, implicit=False)
pi_ = run_standing_wave(om, gamma, dt, 200, implicit=True)
print(f"dt=1.5*dt_sigma  명시적 최대진폭={pe:.1e}  음함수 최대진폭={pi_:.3f}")
# 출력 예: 명시적 최대진폭=3.6e+07  음함수 최대진폭=1.000

코드가 넣는 값은 $\sigma=0.07$ , 밀도 합 $1001.2$ , 길이 $1$ 이다. $1001.2$ 는 물과 공기를 더한 자리고, $0.07$ 은 표면장력이다. $N=64$ 이면 $\Delta x=1/64$ 이고, 같은 식이 주는 $\Delta t_\sigma$ 는 약 $9.32\times 10^{-2}$ 다. $N=128$ 은 $2^{3/2}\approx 2.828$ 로 나눈 약 $3.29\times 10^{-2}$ , $N=256$ 은 다시 나눈 약 $1.16\times 10^{-2}$ 다. 64에서 256으로 간격만 4분의 1이 되는데, 같은 물리 시간의 스텝 수는 8배가 된다. $4^{3/2}=8$ 이라서다. omega*dt_sigma는 해상도와 상관없이 2.221이다. 주석의 예는 $\Delta t=1.5\Delta t_\sigma$ , 200스텝, $\gamma=0.1\omega$ 에서 명시적 최대진폭 $3.6\times 10^{7}$ , 음함수 $1.000$ 이다. 그때 $\omega\Delta t=1.5\times 2.221=3.3315$ 라 심플렉틱 오일러의 2를 넘는다.

$|G|$ 곡선의 가로축이 $\Delta t/\Delta t_\sigma$ 다. mesh와 $\sigma$ 슬라이더가 $\Delta t_\sigma$ 자체를 움직인다.

Δt / Δt_σ = 0.90mesh L/Δx = 128σ = 0.07

explicit |G| = 1.000implicit |G| = 0.447Δt_σ ≈ 3.29e-2 (water/air, L=1)

명시적(주황)은 $\Delta t/\Delta t_\sigma \approx 0.9$ 까지 $|G|=1$ 로 버티다 그 위에서 1을 넘는다. 음함수(시안)는 어떤 스텝에서도 $|G|<1$ 이다. $\sigma$ 를 4배로 올리면 $\Delta t_\sigma$ 는 $\sqrt{1/4}=1/2$ 이 된다. 족쇄의 기울기가 격자에서만 나오는 게 아니다.

완전 결합이 핵심 — 색함수·압력·속도를 한 행렬에#

음함수의 어려움은 표면장력 힘 $\mathbf{f}_\sigma = \sigma\kappa\nabla\psi$ 가 색함수(colour function, 유체 종류를 0~1로 표시하는 VOF 변수) $\psi$ 에 비선형으로 얽힌다는 데 있다. 곡률 $\kappa$ 도 $\psi$ 의 미분으로 정해진다. 논문은 이 항을 새 시각의 $\psi$ 에 대해 Newton 선형화하고, 계면 이류 방정식을 연속·운동량 방정식과 함께 하나의 선형계 $\mathbf{A}\boldsymbol{\phi}=\mathbf{b}$ 로 푼다. 여기서 $\boldsymbol{\phi}=(p, u, v, w, \psi)$ 가 통째로 음함수로 결합된다. Denner 등의 선행 알고리즘은 밀도비 1에 묶여 있었다. Janodet의 기여는 연속·운동량을 보존형으로 쓰고 천이 항에서 밀도를 $\psi$ 에 대해 음함수로 둬, 물/공기 1000:1 같은 큰 밀도비에서도 결합을 유지한 점이다. 계면 이류는 CICSAM 대신 THINC/QQ로 이산화해 큰 CFL에서도 계면을 또렷이 지킨다. OpenFOAM의 분리형 PIMPLE와 달리, 이건 블록 결합(block-coupled) 솔버라야 의미가 산다.

재현하며 부딪힌 것들#

"무조건 안정"이라는 말에는 단서가 붙는다. 음함수 표면장력이 모세관 제약을 깨도, 안정한 최대 $\Delta t$ 는 여전히 유한하다. Galusinski–Vigneaux의 점성-모세관 상한이 그것이다. 밀도비 1000에서는 이 상한이 밀도비 1일 때보다 한 자릿수 더 빡빡해, 논문도 비점성 극한에서 약 $1.5\,\Delta t_\sigma$ 까지만 안정하다고 보고한다. 위의 진동자 모델은 이 점성-모세관 결합을 단순화한 것이라 음함수 쪽이 실제보다 관대하게 보인다는 한계도 있다. 또 후방 오일러는 안정한 대신 빠른 모세관파를 인위적으로 감쇠시킨다. 에너지를 보존해야 하는 문제에서는 이 수치 감쇠가 물리적 감쇠로 오인될 수 있어, 논문은 힘 균형과 에너지 보존을 따로 검증한다.

1.5 Δt_σ에서 끊을지, 밀도비 1000을 버릴지#

음함수가 모세관 CFL을 넘겨도 상한은 남는다. 비점성 극한에서 논문이 말한 안정 구간은 약 $1.5\,\Delta t_\sigma$ 다. 주석의 1.5배, 200스텝에서 명시적 피크가 $3.6\times 10^{7}$ 로 튀고 음함수가 1.000에 머문 것은 진동자 이야기다. 밀도비 1000에서는 Galusinski–Vigneaux 상한이 밀도비 1보다 한 자릿수 더 빡세다. 진동자가 보여주는 관용은 솔버에서 그대로 나오지 않는다.

후방 오일러는 그 대가로 빠른 모세관파를 수치로 죽인다. 에너지를 지켜야 하면 그 감쇠를 물리로 착각하기 쉽다. $\Delta t$ 를 더 민 뒤의 창은 다섯 배에서 열리고 열 배에서 닫힌 기록에 있다.