압력이 격자무늬로 진동할 때 — SIMPLE과 Rhie-Chow 압력-속도 연계

α_p를 1로 두고 압력을 통째로 갈아끼우면 보정 루프가 과보정으로 터진다.

속도장은 멀쩡해 보였다. 압력장을 열었더니 셀마다 위아래로 튀었다. 격자무늬(checkerboard)다. 잔차는 줄어드는데 압력만 톱니다. 버그가 아니었다. 압력과 속도를 같은 격자점에 둔 순간 예고된 결과다.

Naive와 Rhie-Chow 토글을 바꾸면, 심어 둔 톱니 진폭이 그대로 남거나 매끄러운 곡선으로 가라앉는다.

Naive interpolation happened to settle — nudge the grid to see the sawtooth return. (sweeps: 0)

Cells N: 32

Naive에서는 진폭이 안 줄어든다. Rhie-Chow로 바꾸면 같은 격자무늬가 내려앉는다. Reseed checkerboard는 톱니를 다시 심는다.

톱니가 왜 남는지, 연속방정식이 압력 방정식이 되는 자리는 그다음이다. 압축성까지 한 solver로 잇는 자리도 뒤에 있다.

연속방정식은 옆 칸 압력을 보지 못한다#

코로케이트 격자(속도·압력을 같은 셀 중심에 저장)에서 반이산 운동량 방정식은 이렇게 쓴다.

a_P \mathbf{u}_P = \sum_{nb} a_{nb}\mathbf{u}_{nb} - V_P (\nabla p)_P = \mathbf{H}(\mathbf{u}) - V_P (\nabla p)_P

$a_P$ 는 대각 계수, $\mathbf{H}(\mathbf{u})$ 는 이웃 기여, $V_P$ 는 셀 부피, $(\nabla p)_P$ 는 셀 중심 압력구배다.

문제는 셀 중심의 압력구배다. 중심차분으로 쓰면 $(\partial p/\partial x)_P = (p_E - p_W)/2\Delta x$ 가 된다. 여기엔 $p_P$ 가 없다. 셀은 자기 압력을 보지 못하고 두 칸 건너 이웃만 본다. 그래서 홀수 셀과 짝수 셀이 서로 분리된다(odd-even decoupling). $(+,-,+,-,\dots)$ 톱니 압력은 이산 방정식의 영공간(null space)에 숨어 아무 잔차도 만들지 않는다.

Rhie–Chow: 운동량으로 면속도를 만든다#

해법은 면속도를 셀 속도의 단순 평균으로 얻지 않는 것이다. 면에서 운동량 방정식을 다시 적용해, 그 면의 이웃 압력차를 직접 넣는다.

u_f = \overline{\left(\frac{\mathbf{H}(\mathbf{u})}{a_P}\right)}_f - \left(\frac{V}{a_P}\right)_f \frac{p_E - p_P}{\Delta x}

윗줄 막대는 선형보간이다. 뒤쪽 항은 면에서 다시 계산한 압력구배다. 이제 면속도가 $p_E - p_P$ , 즉 바로 옆 셀 압력차에 직접 걸린다. 톱니 모드는 더 이상 보이지 않는 곳에 숨을 수 없다.

면의 $p_E-p_P$ 가 0이 아닌 한 톱니는 그 항에 걸린다. 압력을 운동량에서 지우는 다른 길은 Chorin의 투영법이다.

SIMPLE: 연속방정식을 압력 방정식으로 바꾼다#

압력을 직접 푸는 방정식이 없다는 게 비압축성의 근본 난점이다. SIMPLE(Semi-Implicit Method for Pressure Linked Equations)은 연속방정식 $\nabla\cdot\mathbf{u}=0$ 을 압력 방정식으로 바꾼다. Rhie-Chow 면속도를 연속방정식에 대입하면 압력에 대한 타원형(Poisson) 방정식이 나온다.

\nabla\cdot\left[\left(\frac{V}{a_P}\right)_f \nabla p\right] = \nabla\cdot\left(\frac{\mathbf{H}(\mathbf{u})}{a_P}\right)_f

좌변은 라플라시안(확산꼴), 우변은 예측 속도의 발산이다. 절차는 예측자-수정자다.

추측한 $p^*$ 로 운동량을 풀어 $\mathbf{u}^*$ 예측.
위 압력 방정식을 풀어 $p$ 갱신.
새 압력구배로 면속도·셀속도를 수정해 연속성을 맞춤.
수렴할 때까지 1–3 반복.

PISO는 여기서 수정 단계를 2–3번 더 도는 예측자-수정자로, 과도(transient) 계산에 쓴다. PIMPLE은 큰 시간간격에서 바깥 SIMPLE 루프와 안쪽 PISO 루프를 겹친다.

언더릴랙세이션 없이는 발산한다#

문제가 있다. 스텝마다 압력을 완전히 갱신하면( $p \leftarrow p^*$ ) 커플링이 과보정되어 발산한다. SIMPLE은 압력을 조금씩만 반영한다.

p^{\text{new}} = p^{\text{old}} + \alpha_p\,(p^* - p^{\text{old}})

$\alpha_p$ 는 압력 언더릴랙세이션 계수(0~1)다. 속도에도 $\alpha_u$ 를 따로 건다. 너무 작으면 수렴이 기어가고, 너무 크면 진동하다 터진다. 경험칙은 $\alpha_p \approx 0.3$ , $\alpha_u \approx 0.7$ 이며 둘의 합이 대략 1이 되게 맞춘다.

아래 실습에서 이완 계수를 직접 바꿔보자. 압력-수정 방정식을 이완된 Gauss-Seidel로 푼다.

converging — sweep 0, residual 0.0e+0.

Relaxation factor w (≈ SIMPLE α_p): 1.50

계수를 0.3으로 낮추면 잔차 곡선이 완만하게 기어 내려간다. 1.5 근처에서 가장 빠르다. 2에 붙이면 매 스윕 과보정이 쌓여 잔차가 위로 튄다. SIMPLE의 $\alpha_p$ 도 정확히 이 균형 위에 있다.

하나의 solver로 모든 속도를 — 압축성 확장#

압력-기반 방법의 진짜 매력은 마하수를 안 가린다는 것이다. 압축성에선 연속방정식이 밀도 방정식이 되고, 상태방정식(EOS)으로 밀도를 압력에 건다. $\psi \equiv \partial\rho/\partial p$ (압축률, $=1/c^2$ )를 쓰면 압력 방정식에 시간항이 붙는다.

\frac{\partial (\psi\, p)}{\partial t} + \nabla\cdot(\rho\,\mathbf{u}_f) = 0

이 방정식은 대류-확산 성격을 동시에 가진다. $M\to 0$ 이면 $\psi$ 가 커져 시간항이 지배하지 않고, 압력은 타원형(순간적 전역 커플링)으로 푼다. $M$ 이 크면 쌍곡형(음향파 유한속도 전파)으로 넘어간다. 밀도-기반 solver는 저마하에서 밀도-압력 커플링이 약해 굳지만, 압력-기반은 EOS로 그 커플링을 명시적으로 유지한다. 그래서 저·초음속을 한 코드로 다룬다.

Python으로 체커보드를 되살리고 없앤다#

주장만으로는 부족하다. 두 스텐실을 1차원 주기 격자에서 직접 돌려, 체커보드가 남는지 사라지는지 재본다.

import numpy as np
 
def poisson_sweep(p, f, naive):
    """1D 주기 격자에서 -p'' = f 의 이완 없는 Gauss-Seidel 한 스윕."""
    N = len(p)
    for i in range(N):
        if naive:                       # naive 선형보간: 두 칸 건너뛰는 분리 스텐실
            p[i] = 0.5 * (p[(i - 2) % N] + p[(i + 2) % N] + f[i])
        else:                           # Rhie-Chow: 옆 칸을 잇는 조밀 3점 스텐실
            p[i] = 0.5 * (p[(i - 1) % N] + p[(i + 1) % N] + f[i])
    p -= p.mean()                       # 압력은 상수만큼 자유 → 평균을 고정
    return p
 
def checkerboard_metric(p):
    """(+,-,+,-,...) 성분의 크기. 0이면 톱니 없음."""
    signs = (-1.0) ** np.arange(len(p))
    return abs(np.dot(p, signs)) / len(p)
 
N = 48
x = 2 * np.pi * np.arange(N) / N
f = np.sin(x) + 0.4 * np.sin(2 * x)
f -= f.mean()
 
for naive in (True, False):
    p = 0.8 * (-1.0) ** np.arange(N)    # 체커보드로 초기화
    for _ in range(4000):
        p = poisson_sweep(p, f, naive)
    tag = "naive linear" if naive else "Rhie-Chow  "
    print(f"{tag}  checkerboard = {checkerboard_metric(p):.2e}")
 
# naive linear   checkerboard = 8.00e-01   ← 톱니가 그대로 남는다
# Rhie-Chow      checkerboard = 3.1e-16    ← 기계 정밀도까지 사라진다

같은 소스, 같은 초기 조건, 같은 반복 횟수다. 스텐실만 바꿨다. Naive는 톱니를 못 없애고, Rhie-Chow는 완전히 지운다. 앞 화면과 같은 자리다.

초기 압력은 $0.8\times(-1)^i$ 다. 지표가 $p$ 와 같은 부호열의 내적을 칸수로 나눈 값이면, 처음엔 0.8이다. 48칸, 4000스윕 뒤 naive는 $8.00\times 10^{-1}$ 로 그대로다. 한 점도 안 깎였다. Rhie-Chow는 $3.1\times 10^{-16}$ 이다. 기계정밀도다.

naive가 0.8을 지키는 이유는 스텐실이다. 톱니에서 $i$ 와 $i\pm 2$ 는 부호가 같다. $p_i=0.5(p_{i-2}+p_{i+2}+f_i)$ 는 그 모드를 평균으로 지우지 못한다. 셀 중심 구배 $(p_E-p_W)/2\Delta x$ 도 같다. $E$ 와 $W$ 가 둘 다 $P$ 의 반대 부호라 $p_E-p_W=0$ 이다. 진폭 0.8의 톱니는 중심에서 구배가 0이다. 잔차에 안 잡힌다. 면의 $p_E-p_P$ 는 부호가 반대라 점프 크기가 1.6이다. Rhie-Chow가 잡는 자리다. 소스 $f=\sin x+0.4\sin 2x$ 의 매끈한 부분은 양쪽 스텐실이 다 푼다. 남는 차이는 톱니뿐이다.

α_p를 1로 두면 보정이 먼저 터진다#

$p\leftarrow p^*$ 처럼 압력을 한 스텝에 통째로 갈아끼우면 커플링이 과보정된다. 식은 $\alpha_p$ 만큼만 통과시킨다. 경험칙 $\alpha_p\approx 0.3$ , $\alpha_u\approx 0.7$ 은 합이 1이다. 실습 곡선에서 0.3은 잔차가 기어가고, 1.5 근처가 가장 빠르며, 2에 붙이면 스윕마다 잔차가 다시 튄다. SIMPLE의 $\alpha_p$ 는 0과 1 사이에 두고, 그 곡선의 왼쪽 편에 앉힌다.

이완을 바꿨더니 수렴해가 달라지는 자리는 Rhie–Chow가 언더릴랙세이션을 기억하는 면속도다. 코로케이트에서 면속도를 셀 속도의 평균으로 두면 옆 칸 압력차가 빠지고, 톱니는 잔차 밖에 남는다. $\psi=\partial\rho/\partial p$ 로 시간항을 살리는 일은 그다음이다. 톱니가 살아있는 압력 방정식에 압축성을 얹으면 음속만 빨라지고 격자무늬는 그대로다.