압력을 방정식에서 지운 사람 — Chorin의 투영법과 분수 시간전진

유입과 유출이 한 방울만 어긋나도 Neumann Poisson은 해가 없다. 투영을 끄면 그 발산이 스텝마다 남는다.

1967년, Alexandre Chorin은 압력을 방정식에서 잠시 지워버렸다. 비압축성 Navier–Stokes를 풀 때 가장 다루기 까다로운 항이 압력이었기 때문이다. 압력에는 시간 미분이 없다. 밀도가 상수라 상태방정식으로 압력을 구할 수도 없다. Chorin의 답은 대담했다. 일단 압력을 무시하고 속도를 전진시킨 뒤, 그 결과를 강제로 비발산(divergence-free) 상태로 되돌리는 것이다. 이 글은 그 분수 시간전진(fractional step)—흔히 투영법(projection method)—이 어떻게 작동하는지, Helmholtz–Hodge 분해라는 한 줄의 정리에서 출발해 직접 굴러가는 2D 솔버까지 만들어 본다. 끝까지 읽으면 왜 비압축성 코드의 실행 시간 대부분이 압력 Poisson 방정식에 잡아먹히는지 알게 된다.

왼쪽 패널의 소스 세기를 0에서 14로 올리면 발산 색이 살아나고, 오른쪽 투영은 소용돌이만 남긴다.

source strength6.0max|∇·u*| = 0.000 → 0.0000

Red = positive divergence (source), blue = negative (sink), dark = zero. Raise the source and the left panel lights up; the right panel stays dark — projection strips the compressible part and keeps only the swirl.

소스가 0이면 양쪽이 거의 같다. 14까지 올리면 왼쪽만 빨강·파랑으로 물들고, 오른쪽은 어둡게 남는다. 투영이 압축 성분을 통째로 걷어낸 결과다.

압력에는 시간 미분이 없다#

비압축성 Navier–Stokes를 적으면 문제가 눈에 보인다.

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}, \qquad \nabla\cdot\mathbf{u} = 0

여기서 $\mathbf{u}$ 는 속도, $p$ 는 압력, $\nu$ 는 동점성계수(운동량 확산율)다. 운동량 방정식은 $\mathbf{u}$ 의 시간 발전을 알려준다. 그런데 두 번째 식 $\nabla\cdot\mathbf{u}=0$ 은 발전 방정식이 아니다. 매 순간 지켜야 할 구속 조건이다.

압축성 유동에서는 연속방정식이 밀도의 발전을 주고, 밀도가 상태방정식을 통해 압력을 정한다. 비압축성에서는 그 사슬이 끊긴다. 압력은 시간에 따라 "발전"하지 않는다. 압력은 속도가 발산하지 않도록 매 순간 스스로를 맞추는 라그랑주 승수(constraint를 강제하는 미지수)일 뿐이다. 그러니 압력을 시간 전진시키려는 시도 자체가 헛일이다.

모든 속도장은 두 조각이다 — Helmholtz–Hodge#

돌파구는 벡터 미적분의 오래된 정리에 있다. 적당한 영역에서 임의의 벡터장 $\mathbf{w}$ 는 비발산 부분과 기울기 부분의 합으로 유일하게 쪼개진다.

\mathbf{w} = \mathbf{u} + \nabla\phi, \qquad \nabla\cdot\mathbf{u} = 0

$\mathbf{u}$ 는 소용돌이(비발산) 성분, $\nabla\phi$ 는 압축(기울기) 성분이다. 이 분해에서 $\mathbf{u}$ 만 뽑아내는 연산을 투영 연산자 $\mathbb{P}$ 라 부른다. 방법은 간단하다. 위 식의 발산을 취하면 $\nabla\cdot\mathbf{u}=0$ 이므로

\nabla^2\phi = \nabla\cdot\mathbf{w}

이 Poisson 방정식으로 $\phi$ 를 풀고, 원래 장에서 $\nabla\phi$ 를 빼면 비발산 부분만 남는다.

\mathbf{u} = \mathbb{P}\mathbf{w} = \mathbf{w} - \nabla\phi

예측하고, 투영한다 — 분수 시간전진#

Chorin의 알고리즘은 이 분해를 시간 적분에 그대로 얹는다. 한 스텝을 두 조각으로 나눈다.

예측자(predictor) — 압력을 뺀 채 속도를 전진시켜 중간 속도 $\mathbf{u}^*$ 를 얻는다.

\frac{\mathbf{u}^* - \mathbf{u}^n}{\Delta t} = -(\mathbf{u}^n\cdot\nabla)\mathbf{u}^n + \nu\nabla^2\mathbf{u}^n

$\mathbf{u}^*$ 는 비발산이 아니다. 압력을 무시했으니 당연하다.

투영자(corrector) — $\mathbf{u}^*$ 를 비발산으로 되돌린다. 압력 Poisson을 풀고,

\nabla^2\phi = \frac{1}{\Delta t}\nabla\cdot\mathbf{u}^*

기울기를 빼서 다음 스텝 속도를 만든다.

\mathbf{u}^{n+1} = \mathbf{u}^* - \Delta t\,\nabla\phi

$\phi$ 가 사실상 압력 역할을 한다( $p^{n+1}\approx\rho\phi$ ). 압력은 발전시키는 대상이 아니라, 매 스텝 구속 조건을 맞추려고 새로 푸는 값이다. 이 두 줄이 투영법의 전부다.

압력 Poisson과 경계조건의 함정#

여기서 실무의 고통이 시작된다. 예측자는 값싸다. 그러나 투영자는 매 스텝 Poisson 방정식을 하나 푼다. 큰 격자에서 이 타원형 solve가 전체 비용의 70~90%를 차지한다. GMRES나 멀티그리드가 비압축성 코드의 심장인 이유다.

경계조건도 함정투성이다. $\phi$ 의 경계조건은 벽에서 Neumann( $\partial\phi/\partial n = 0$ )이다. Neumann 문제는 상수만큼의 불확정성이 있고, 해가 존재하려면 적합성 조건—영역 전체에서 $\int \nabla\cdot\mathbf{u}^*\,dV = 0$ —을 만족해야 한다. 유입·유출 유량이 안 맞으면 Poisson 방정식은 애초에 풀리지 않는다. NaN이 뜨면 경계 유량부터 의심하라.

collocated(비엇갈림) 격자에서는 압력과 속도를 같은 점에 두면 체커보드 진동이 생긴다. 엇갈림 격자를 쓰거나 Rhie–Chow 보간으로 막아야 한다. 이완 계수를 그 면속도가 기억하면, 수렴해 자체가 달라진다. 참고로 SIMPLE 계열은 같은 투영 아이디어를 매 스텝 여러 번 반복해 정상상태로 몰고 가는 반면, 여기 분수 시간전진은 스텝당 한 번의 투영으로 비정상 유동을 시간 정확하게 따라간다.

Python: 이중 전단층을 FFT로 굴린다#

주기 경계에서는 압력 Poisson을 FFT로 단번에 푼다. Fourier 공간에서 라플라시안이 $-|\mathbf{k}|^2$ 곱으로 바뀌기 때문이다. 이중 전단층(두 개의 맞물린 제트)을 초기조건으로 주면 Kelvin–Helmholtz 소용돌이로 말려 올라간다. 투영법의 고전 검증 문제다.

import numpy as np
 
N = 128                # 격자 (한 변)
h = 1.0 / N            # 격자 간격
nu = 5e-3              # 동점성 → Re = U L / nu ≈ 200
dt = 2e-3              # 확산·이류 안정 범위 안
steps = 3000
 
x = (np.arange(N) + 0.5) * h
X, Y = np.meshgrid(x, x, indexing='ij')
 
# 이중 전단층 + 작은 교란
rho, delta = 1.0 / 30.0, 0.05
u = np.where(Y <= 0.5, np.tanh((Y - 0.25) / rho), np.tanh((0.75 - Y) / rho))
v = delta * np.sin(2 * np.pi * X)
 
# FFT Poisson용 파수
k = 2 * np.pi * np.fft.fftfreq(N, d=h)
KX, KY = np.meshgrid(k, k, indexing='ij')
K2 = KX**2 + KY**2
K2[0, 0] = 1.0         # 평균 모드 0으로 나눔 방지
 
def divergence(a, b):
    dadx = (np.roll(a, -1, 0) - np.roll(a, 1, 0)) / (2 * h)
    dbdy = (np.roll(b, -1, 1) - np.roll(b, 1, 1)) / (2 * h)
    return dadx + dbdy
 
def projection_correct(a, b):
    # Laplacian(phi) = div  →  a <- a - grad(phi)
    phi_hat = np.fft.fft2(divergence(a, b)) / (-K2)
    phi_hat[0, 0] = 0.0
    phi = np.real(np.fft.ifft2(phi_hat))
    dpx = (np.roll(phi, -1, 0) - np.roll(phi, 1, 0)) / (2 * h)
    dpy = (np.roll(phi, -1, 1) - np.roll(phi, 1, 1)) / (2 * h)
    return a - dpx, b - dpy
 
def predictor(a, b):
    # 이류(중심차분) + 확산, 명시적 Euler 한 스텝
    ax = (np.roll(a, -1, 0) - np.roll(a, 1, 0)) / (2 * h)
    ay = (np.roll(a, -1, 1) - np.roll(a, 1, 1)) / (2 * h)
    bx = (np.roll(b, -1, 0) - np.roll(b, 1, 0)) / (2 * h)
    by = (np.roll(b, -1, 1) - np.roll(b, 1, 1)) / (2 * h)
    lap = lambda f: (np.roll(f, -1, 0) + np.roll(f, 1, 0)
                     + np.roll(f, -1, 1) + np.roll(f, 1, 1) - 4 * f) / h**2
    astar = a + dt * (-(a * ax + b * ay) + nu * lap(a))
    bstar = b + dt * (-(a * bx + b * by) + nu * lap(b))
    return astar, bstar
 
for n in range(steps):
    us, vs = predictor(u, v)                 # 예측: 중간 속도 u*
    d_before = np.abs(divergence(us, vs)).max()
    u, v = projection_correct(us, vs)        # 투영: 발산 제거
    if n % 500 == 0:
        d_after = np.abs(divergence(u, v)).max()
        print(f"step {n:4d}  |div u*|={d_before:.2e} -> |div u|={d_after:.2e}")

출력을 보면 매 스텝 $|\nabla\cdot\mathbf{u}^*|$ 는 $10^{-2}$ 규모인데 투영 뒤 $|\nabla\cdot\mathbf{u}|$ 는 $10^{-13}$ 까지 떨어진다. FFT 투영이 발산을 기계 정밀도로 죽인 것이다.

격자 한 변이 128이면 간격은 1/128이다. 동점성 5e-3, 속도 규모 1, 영역 길이 1이면 Re는 1을 5e-3으로 나눈 200이다. 주석의 Re 약 200과 같다. 시간 간격 2e-3를 3000번 쌓으면 적분 시간은 6이다. 전단층 두께 1/30은 약 0.033, 간격 1/128이 약 0.0078이니 네 칸 남짓이다. 교란 진폭 0.05는 제트 속도 1의 5%라서, 말림은 이 작은 사인파에서 시작된다. 예측 직후 발산 $10^{-2}$ 와 투영 뒤 $10^{-13}$ 의 비는 열한 자릿수다. FFT가 평균 모드를 0으로 고정하는 한 줄은 Neumann 문제의 상수 자유도를 닫는다. 유입과 유출이 안 맞아 적합성 조건이 깨지면, 그 분모만 손봐서는 해가 생기지 않는다. 비용 70%에서 90%라는 말도 같은 행렬을 가리킨다. 한 스텝에서 예측자가 가져가는 몫은 나머지 10%에서 30%뿐이다.

시간 간격 2e-3에서 보정은 중간 속도에서 그 간격만큼 $\phi$ 의 기울기를 빼는 일이다. 발산 $10^{-2}$ 를 지우려면 $\phi$ 의 라플라시안은 발산을 시간 간격으로 나눈 값, 대략 5가 된다. 압력 구배가 한 스텝에 속도를 그만큼 되돌린다. 투영을 끄면 이 보정이 속도에 들어가지 않고 발산으로 남는다. 3000스텝은 그 찌꺼기가 격자를 채우기에 넉넉하다.

투영을 끄면 벌어지는 일#

투영자가 없으면 속도장은 매 스텝 조금씩 발산을 쌓는다. 질량이 보존되지 않고, 소용돌이는 뭉개지며, 곧 격자를 채운 잡음이 된다. 아래 시뮬레이션에서 직접 확인해보자.

dt0.60max|∇·u| = 0.000

Double shear layer rolling up. Red/blue = vorticity sign. Turn projection OFF and max|∇·u| climbs while the vortices dissolve into noise.

Projection을 ON으로 두면 전단층이 깔끔한 소용돌이 쌍으로 말려 올라간다. OFF로 바꾸는 순간 우측 상단의 max|∇·u|가 치솟고 소용돌이 무늬가 무너진다. dt를 키우면 이류가 더 공격적이 되어 무너지는 속도가 빨라진다. 이 한 번의 토글이 "왜 투영이 필요한가"에 대한 가장 짧은 대답이다.

유량이 한 방울만 어긋나면 Poisson은 답이 없다#

비압축에서 압력은 발전시키는 변수가 아니라, 매 스텝 $\nabla\cdot\mathbf{u}=0$ 을 강제하는 라그랑주 승수다. 투영은 압력을 뺀 채 $\mathbf{u}^*$ 를 만든 뒤 Poisson으로 $\phi$ 를 풀어 기울기를 뺀다. 그 순서 자체는 짧다. 깨지는 자리는 경계다. 유입과 유출이 안 맞으면 적합성 조건이 없고, 행렬은 상수 자유도만 남긴 채 해를 내놓지 않는다.

한 스텝의 투영으로 시간 정확도를 지키려면 대류를 그 스텝 안에서 다시 풀지 않는다. 점성만 차수를 올리고 대류는 바깥으로 미는 선택은 Gear가 대류를 외삽하다 오차를 키운 자리에서 갈린다. 예측자의 대류를 더 정확해 보이게 만들다 보면, 투영이 맞춰 둔 발산 제로가 흔들린다.