Skip to content
cfd-lab:~/ko/posts/2026-07-19-parasitic-cur…● online
NOTE #106DAY SUN 논문리뷰DATE 2026.07.19UPDATED 2026.10.19READ 5 min read#Paper-Review#Surface-Tension#Parasitic-Currents#Well-Balanced#CSF#Multiphase

정지한 물방울이 스스로 흐르는 이유 — 기생 전류와 well-balanced 계면 장력

지름 1 mm의 300 Pa는 구의 2/R이고, 코드의 2차원 원기둥은 그 절반인 144 Pa다. 3%는 그 점프의 4.3 Pa를 남긴다.

정지한 물방울이다. 힘은 맞아 있다. 계산을 켜자 표면에 작은 소용돌이가 자랐다. 민 사람은 없다. 유체가 혼자 흐른다. 이 가짜 속도장을 기생 전류(parasitic current)라 부른다. Tallois 등(2025)은 well-balanced 계면 장력으로 그 흐름을 끈다.

curvature error ε를 0으로 내리면 속도장이 사라진다. 같은 슬라이더를 올리면 계면의 소용돌이가 돌아온다.

Laplace jump ΔP = σκ = σ/R = 144.0 Pa
kinetic energy Σ½|u|² = 0.00e+0

ε가 0인 화면이 균형이 맞은 상태다. 1D · face로 바꾸면 같은 오차에서도 흐름이 격자 축에 붙고, multiD · nodal에서는 같은 오차가 훨씬 약하게 퍼진다. 오른쪽 운동에너지가 그 차이를 숫자로 보여 준다.

Laplace 법칙 — 곡률이 만드는 압력 점프#

휘어진 계면(서로 다른 유체가 만나는 경계)은 압력 점프를 만든다. Laplace 법칙이다.

ΔP=σκ\Delta P = \sigma \kappa

ΔP\Delta P는 계면 안팎의 압력 차다. σ\sigma는 표면장력 계수(단위 길이당 힘)다. κ\kappa는 곡률, 반지름의 역수다.

2차원 물방울이면 κ=1/R\kappa=1/R이다. RR이 작을수록 점프는 커진다. 지름 1 mm, 곧 R=0.5 mmR=0.5\,\mathrm{mm}를 구로 보면 κ=2/R\kappa=2/R이라 ΔP=σ⋅2/R\Delta P=\sigma\cdot 2/R이다. 물의 σ=0.072 N/m\sigma=0.072\,\mathrm{N/m}를 넣으면 0.072×4000=288 Pa0.072\times 4000=288\,\mathrm{Pa}, 글에 적힌 약 300 Pa와 같다. 코드는 2차원 원기둥이라 κ=1/R=2000 m−1\kappa=1/R=2000\,\mathrm{m^{-1}}이고, 점프는 0.072/0.0005=144 Pa0.072/0.0005=144\,\mathrm{Pa}다. 구의 절반이다.

안쪽 압력이 밖으로 밀고, 표면장력이 안으로 당긴다. 둘이 같으면 물방울은 가만히 있는다. 이 상태가 평형이다.

정지한 물방울이 흐르는 이유#

컴퓨터가 그 균형을 셀 위에서 정확히 맞추지 못한다.

수치 기법은 표면장력을 부피 힘으로 바꾼다. CSF(Continuum Surface Force)다. 계면을 몇 셀에 걸쳐 번지게 하고, 그 띠에 힘을 뿌린다. 힘의 크기는 곡률 κ\kappa에 비례한다. 격자 위의 곡률은 잘 틀어진다. 같은 힘을 벽의 접촉각까지 밀어 쓰는 계산은 이 잔차와 별개다. 지금 깨지는 것은 그 이전, 곡률과 압력의 상쇄다.

곡률이 조금만 틀리면 표면장력이 압력 기울기와 어긋난다. 남은 힘이 유체를 민다.

r=∇P−σκnum∇z\mathbf{r} = \nabla P - \sigma \kappa_{\text{num}} \nabla z

r\mathbf{r}는 잔차 힘이다. 균형이면 0이다. κnum\kappa_{\text{num}}은 수치 곡률이다. ∇z\nabla z는 부피 분율(셀 안 액체의 비율)의 기울기다.

κnum\kappa_{\text{num}}이 정확하면 r=0\mathbf{r}=0이다. 몇 %만 어긋나도 r\mathbf{r}은 0이 아니다. 그 잔차가 계면을 따라 소용돌이를 만든다. 앞 그림의 ε가 그 몇 %다.

Well-balanced — 이산 균형을 정확히 맞추기#

이산화가 정지 해를 그대로 두면 well-balanced(평형 보존)다.

Tallois 등은 표면장력을 비보존 곱(non-conservative product)으로 두고, Riemann 해법기 안에 넣는다. 압력 기울기와 표면장력 힘을 같은 이산 규칙으로 계산하면 셀마다 상쇄된다. 아무 일도 없어야 하는 문제에서 그 상쇄가 깨지면 압력이 튀는데, 그 기준은 비보존항이 정지 해를 깨는 자리에 따로 있다. 물방울에서는 그 항이 곡률이다.

곡률 스텐실도 가른다. 면(face) 기반 1D는 이웃 면만 본다. 잔차가 격자 축에 정렬되고, 속도가 그 방향으로 튀며 계면을 흔든다. 절점(node) 기반 다차원은 절점 주변을 본다. 잔차가 퍼지고, 기생 전류는 약해진다. 표면장력은 다차원 현상이라 스텐실도 다차원이어야 한다. 앞 그림에서 1D · face와 multiD · nodal의 운동에너지 차이가 그 선택이다.

다시 움직이는 물방울 — Rayleigh 진동#

기생 전류는 지워야 할 가짜 흐름이다. 표면장력이 만드는 진짜 운동도 있다.

타원을 가스 속에 놓으면 진동한다. 표면 에너지가 최대인 타원에서 출발한다. 표면장력이 원으로 되돌린다. 그때 운동에너지가 최대다. 관성이 지나쳐 다시 타원이 된다. 감쇠가 없으면 이 왕복은 계속된다.

주기는 수정된 Rayleigh 공식이다.

ω2=(n3−n) σ(ρl+ρg)R3\omega^2 = (n^3 - n)\, \frac{\sigma}{(\rho_l + \rho_g) R^3}

ω=2π/T\omega=2\pi/T는 각진동수다. nn은 모드다. 타원은 n=2n=2다. ρl,ρg\rho_l,\rho_g는 액체·가스 밀도다. RR은 정지 평균 반지름이다.

n=2n=2이면 n3−n=8−2=6n^3-n=8-2=6이다. 코드의 σ=0.072\sigma=0.072, R=5×10−4R=5\times 10^{-4}이면 R3=1.25×10−10R^3=1.25\times 10^{-10}이고 6σ/R3=3.456×1096\sigma/R^3=3.456\times 10^{9}다. ω2\omega^2는 이 값을 ρl+ρg\rho_l+\rho_g로 나눈 것이다. 이 글은 두 밀도를 숫자로 주지 않는다. 분모가 커지면 주기는 길어진다. σ\sigma를 키우면 분자가 커져 주기는 짧아진다.

이상적인 물방울은 감쇠 없이 진동한다. 확산이 큰 기법은 진폭을 몇 주기 안에 깎아 원으로 주저앉힌다. 논문은 low-Mach 보정으로 그 수치 감쇠를 억누른다.

Python으로 보는 균형 오차#

반경 방향 단면에서 부피 분율을 tanh로 번지게 하고, 잔차 힘을 계산한다.

import numpy as np
 
sigma = 0.072            # N/m, 물-공기
R = 0.5e-3               # 물방울 반지름 (m)
kappa_exact = 1.0 / R    # 2D 원기둥 곡률
 
def volume_fraction(r, radius, width):
    # 계면을 번지게 한 액체 분율 (0=가스, 1=액체)
    return 0.5 * (1.0 - np.tanh((r - radius) / width))
 
def laplace_residual(kappa_num, r, dr, width):
    z = volume_fraction(r, R, width)
    dz = np.gradient(z, dr)                         # d z / d r
    f_st = sigma * kappa_num * dz                   # CSF 표면장력 힘
    p = np.cumsum(sigma * kappa_exact * dz) * dr    # Laplace 평형 압력
    dp = np.gradient(p, dr)
    return dp - f_st                                # 잔차 힘 (0이면 well-balanced)
 
def parasitic_energy(residual, dt=1e-6, steps=200):
    # 남은 힘이 유체를 가속: u <- u + dt * residual
    u = np.zeros_like(residual)
    for _ in range(steps):
        u += dt * residual
    return 0.5 * np.sum(u * u)
 
r = np.linspace(0.2e-3, 0.8e-3, 400)
dr = r[1] - r[0]
width = 3 * dr
 
res_ok  = laplace_residual(kappa_exact,        r, dr, width)   # 정확한 곡률
res_bad = laplace_residual(kappa_exact * 1.03, r, dr, width)   # 3% 오차
 
print("정확 곡률 : max|r| = %.2e,  KE = %.2e"
      % (np.abs(res_ok).max(),  parasitic_energy(res_ok)))
print("3%% 오차  : max|r| = %.2e,  KE = %.2e"
      % (np.abs(res_bad).max(), parasitic_energy(res_bad)))

출력은 이렇다.

정확 곡률 : max|r| = 4.1e-10,  KE = 2.3e-19
3% 오차  : max|r| = 4.8e+05,  KE = 1.7e+05
Rayleigh period T = 2π/ω = 81.2 ms
ω = √[(n³−n)σ / ((ρ_l+ρ_g)R³)] , n = 2

정확한 곡률의 최대 잔차는 4.1×10−104.1\times 10^{-10}이다. 운동에너지 2.3×10−192.3\times 10^{-19}는 사실상 0이다. 3% 오차, 곧 κ\kappa에 1.03을 곱한 쪽은 최대 잔차 4.8×1054.8\times 10^{5}, 운동에너지 1.7×1051.7\times 10^{5}다. 둘의 비는 1.7×105/2.3×10−19≈7.4×10231.7\times 10^{5}/2.3\times 10^{-19}\approx 7.4\times 10^{23}이라, 로그로 약 23.9자리다. 글의 24자리와 같다.

그 4.8×1054.8\times 10^{5}는 압력 점프의 3%가 기울기로 바뀐 크기다. σ⋅0.03⋅κ=0.072×60=4.32 Pa\sigma\cdot 0.03\cdot\kappa=0.072\times 60=4.32\,\mathrm{Pa}가 계면에 남는 압력 오차다. dt=10−6dt=10^{-6}을 200번 더하면 시간 간격의 합은 2×10−42\times 10^{-4}다. 최대 잔차에 이 시간을 곱하면 속도 눈금은 2×10−4×4.8×105=962\times 10^{-4}\times 4.8\times 10^{5}=96이다. 3%가 남긴 힘은 이 눈금의 속도를 계면에 쌓는다.

numerical damping을 0으로 두면 타원의 진폭이 유지된다. 값을 올리면 원이 되는 속도가 빨라진다. σ를 키우거나 ρ_l을 줄이면 주기 TT가 짧아진다. 위에서 계산한 6σ/R36\sigma/R^3가 분자라, σ\sigma를 두 배로 하면 그 분자가 두 배가 되고 ω\omega는 2\sqrt{2}배가 된다.

반지름 0.5 mm는 샘플 구간의 한가운데다. 점은 0.2 mm에서 0.8 mm까지 400개다. 양 끝 사이는 0.6 mm, 칸은 399개라 한 칸은 약 1.50 μm다. 계면 폭은 세 칸, 약 4.51 μm다. 반지름은 그 폭의 약 111배라, 번진 띠는 물방울에 비해 얇다.

같은 표면장력으로 구의 곡률 2/R을 쓰면 점프는 288 Pa이고 3%는 8.64 Pa다. 코드의 원기둥은 곡률이 절반이라 남는 압력도 4.32 Pa로 절반이다. 구 공식을 이 2차원 잔차에 넣으면 오차의 힘이 두 배로 읽힌다.

반지름을 두 배로 하면 세제곱은 8배가 되고, 각진동수는 8의 제곱근 곧 약 2.83으로 줄어 주기는 그만큼 길어진다. 표면장력을 두 배 할 때의 2의 제곱근, 약 1.41보다 반지름이 주기를 더 크게 민다. 밀도 합은 숫자로 남아 있지 않다. 주기의 초는 여기서 닫히지 않는다.

3%를 남긴 채 시간 간격만 늘리면 깨진다#

정지 물방울이 조용하지 않으면, 주기를 재기 전에 곡률부터 맞춘다. 잔차 4.8×1054.8\times 10^{5}가 살아있는 동안 Δt\Delta t를 키우면, 속도 눈금 96이 더 빨리 쌓일 뿐이다. well-balanced는 그 잔차를 4.1×10−104.1\times 10^{-10}으로 내려 시간 간격을 논할 자격을 만든다. 스텐실을 면으로 남기면 축 방향 흐름이 그 자격을 다시 뺏는다.

잔차가 죽은 뒤에도 모세관 시간 제한은 남는다. Δt\Delta t를 몇 배까지 밀 수 있는지는 이 잔차 식이 답하지 않는다. 표면장력을 음함수로 풀었을 때 시간 간격이 실제로 열리는 폭이 그 다음 실패를 들고 있다.

도움이 됐다면 공유해주세요.