[논문 리뷰] 언더릴랙세이션을 없애는 법 — 압력기반 완전연성 솔버의 선형화

고정계수 반복이 발산 대신 다른 근에 앉는 경계는 기울기 절댓값 M/Co=1이다.

수렴하지 않는 계산 앞에서 가장 먼저 손이 가는 손잡이는 언더릴랙세이션 계수다. 0.7에서 0.5로, 다시 0.3으로 내린다. 반복 횟수는 늘고 계산 시간은 두 배가 되지만 어쨌든 답은 나온다.

Denner(2018)는 그 손잡이를 아예 떼어내자고 말한다. 대신 손대야 할 곳은 비선형 항을 자르는 방식, 곧 선형화다. 그 주장을 잔차 숫자로 따라가 봤다. 고정계수 선형화와 Newton 선형화가 이산화된 연속방정식에서 각각 무엇을 남기고 무엇을 지우는지 보고, 지워진 항이 언제부터 치명적이 되는지 장난감 모형으로 눈금을 매긴다.

논문: F. Denner, Fully-coupled pressure-based algorithm for compressible flows: linearisation and iterative solution strategies, arXiv:1807.04232 (2018). Imperial College London. 전마하 압축성 유동을 위한 완전연성(fully-coupled) 압력기반 알고리즘에서, 선형화 전략과 반복 해법이 성능·안정성에 미치는 영향을 체계적으로 비교한 논문.

마하 슬라이더를 0.001에서 3으로 올리면 오른쪽 밀도 쪽 추가 무거워지고, fixed-coefficient 버튼은 그 추를 통째로 지운다.

linearisation

Mach M0.050Courant Co0.400

아래 막대의 흰 표시는 타원형에서 쌍곡형으로 미끄러진다. 저마하에서는 압력-속도 결합이 빔을 붙들고, 마하가 오르면 압력-밀도 결합이 기울기를 가져간다.

초음속에서 멈춘 것은 솔버가 아니라 반복이었다#

압축성 솔버는 크게 두 갈래다. 밀도기반(density-based)은 연속방정식을 밀도의 수송방정식으로 본다. 충격파가 있는 초음속에서 잘 돈다. 대신 저마하에서 무너진다. 마하수가 0으로 가면 밀도와 압력의 결합이 사라지기 때문이다.

압력기반(pressure-based)은 반대다. 연속방정식을 압력 방정식으로 쓰고, 밀도는 상태방정식으로 따로 구한다. 저마하에서 강하다.

문제는 그 사이다. 천음속(transonic·마하수 1 부근) 영역에서는 압력-속도 결합과 압력-밀도 결합이 동시에 세다. 두 비선형성이 겹치는 이 구간에서 압력기반 알고리즘의 수렴이 흔들린다. SIMPLE 같은 분리형(segregated) 해법이 언더릴랙세이션 없이 못 도는 이유이기도 하다.

완전연성 해법은 연속·운동량·에너지 방정식을 하나의 선형계에 넣고 동시에 푼다. 메모리는 더 먹지만 결합이 강해진다. 그런데 "하나의 선형계에 넣는다"는 말 자체가 이미 선택을 포함한다. 원래 방정식은 비선형이다. 무엇을 미지수로 두고 무엇을 계수로 넘길지 정해야 선형계가 만들어진다.

압력은 두 개의 일을 동시에 한다#

압력기반 알고리즘의 성공은 압력이 두 가지 역할을 겸한다는 데서 나온다.

저마하에서 압력은 속도장에 대한 구속조건이다. 연속방정식은 타원형(elliptic·해가 전 영역에 즉시 영향을 주는 성질) 압력 방정식이 된다. 밀도는 거의 상수다.

초음속에서는 반대다. 압력은 밀도와 직접 묶이고, 연속방정식은 쌍곡형(hyperbolic·정보가 유한 속도로 전파되는 성질)에 가까워진다. 압력-속도 결합은 뒷전이 된다.

면(face)에서의 질량유속 $\tilde\rho_f\vartheta_f$ 를 압력으로 미분하면 두 결합의 무게가 그대로 나온다.

\frac{\partial(\tilde\rho_f\vartheta_f)}{\partial p} = \underbrace{\rho\,\hat d}_{\text{속도 쪽}} + \underbrace{\vartheta\,\frac{\partial\rho}{\partial p}}_{\text{밀도 쪽}}

$\vartheta_f$ 는 이송 속도(advecting velocity·면을 통과하는 유속으로, 운동량가중보간으로 얻는다), $\hat d$ 는 그 보간에 들어가는 압력 감쇠 계수다. 이상기체 등온 근사에서 $\partial\rho/\partial p = 1/a_T^2$ 이므로 두 항의 비는 다음 한 숫자로 줄어든다.

\frac{\text{밀도 쪽}}{\text{속도 쪽}} = \frac{\vartheta/a_T^2}{\rho\,\hat d} = \frac{M}{Co}

$M = u_0/a_T$ 는 마하수, $Co = a_T\Delta t/\Delta x$ 는 음향 Courant 수(음파가 한 스텝에 격자 몇 칸을 지나는지)다. 이 비가 1을 넘으면 밀도 쪽이 속도 쪽보다 무겁다. 앞 시소에서 fixed-coefficient가 오른쪽 추를 지운 것은, 그 항을 행렬에 넣지 않는다는 뜻이다.

비선형 항을 자르는 두 가지 방법#

일반적인 비선형 항 $\alpha^{(n+1)}\varphi^{(n+1)}$ 을 생각한다. $n$ 은 비선형 반복 횟수다. 고정계수 선형화(fixed-coefficient·계수를 지연시킨다)는 주 변수만 음함수로 둔다.

\alpha^{(n+1)}\varphi^{(n+1)} \approx \alpha^{(n)}\varphi^{(n+1)}

구현이 쉽다. 계수를 이전 반복 값으로 채우기만 하면 된다.

Newton 선형화는 두 변수 모두에 대해 1차 전개한다.

\alpha^{(n+1)}\varphi^{(n+1)} \approx \alpha^{(n)}\varphi^{(n+1)} + \alpha^{(n+1)}\varphi^{(n)} - \alpha^{(n)}\varphi^{(n)}

항이 하나 늘고, 대신 $\alpha$ 도 음함수로 다뤄야 한다. 압력기반 알고리즘에서 $\alpha$ 가 밀도라면 이는 곧 $\rho = \rho(p,T)$ 를 통해 압력에 대한 음함수 의존성을 행렬에 넣는다는 뜻이다. 압력은 이미 모든 방정식의 주 미지수이므로 새로운 0이 아닌 행렬 성분은 생기지 않는다. 공짜에 가깝다는 것이 이 논문의 중요한 실무적 관찰이다.

연속방정식에서 무엇이 사라지는가#

이산화된 연속방정식에 두 선형화를 각각 적용한다. 고정계수는 이렇게 된다.

\frac{\rho_P^{(n+1)} - \rho_P^{(t-\Delta t)}}{\Delta t}V_P + \sum_f \tilde\rho_f^{(n)}\vartheta_f^{(n+1)}A_f = 0

Newton은 항이 늘어난다.

\frac{\rho_P^{(n+1)} - \rho_P^{(t-\Delta t)}}{\Delta t}V_P + \sum_f \Big(\tilde\rho_f^{(n)}\vartheta_f^{(n+1)} + \tilde\rho_f^{(n+1)}\vartheta_f^{(n)} - \tilde\rho_f^{(n)}\vartheta_f^{(n)}\Big)A_f = 0

$V_P$ 는 셀 부피, $A_f$ 는 면적, 위첨자 $(t-\Delta t)$ 는 이전 시간층이다. 앞 절의 시소가 그대로 여기 들어 있다. $\tilde\rho_f^{(n)}\vartheta_f^{(n+1)}$ 은 속도 쪽, $\tilde\rho_f^{(n+1)}\vartheta_f^{(n)}$ 은 밀도 쪽이다.

고정계수는 밀도 쪽 항을 통째로 버린다. 저마하에서는 문제가 없다. 버린 쪽이 가벼웠기 때문이다. 마하수가 커지면 버린 쪽이 무거운 쪽이 된다. 논문의 표현대로, 고정계수 선형화는 비압축성 압력기반 틀에서 유도된 것이라 큰 마하수에서 성능과 안정성이 극히 제한된다.

음향 항을 음해로 넣어도 시간 간격이 다시 잘게 쪼개지는 수가 있다. 상대속도에 두 번째 $1/M^{2}$ 가 숨어 $\Delta t$ 가 100배 줄어든 전마하 이상유동은, 압력-밀도 결합을 행렬에 넣었다고 강성이 끝나는 것이 아님을 보여 준다.

운동량·에너지 방정식의 네 갈래#

운동량·에너지 방정식의 이송항 $\tilde\rho_f\vartheta_f\tilde\varphi_f$ 에는 미지수가 셋이라 선택지가 넷으로 늘어난다. $\varphi$ 는 속도 성분 또는 비전엔탈피다.

이름	음함수로 두는 것	성격
고정계수	$\tilde\varphi_f$ 만	비압축성 관행 그대로
$\rho$ -Newton	$\tilde\varphi_f$ , $\tilde\rho_f$	압력-밀도 결합을 음함수로
$\vartheta$ -Newton	$\tilde\varphi_f$ , $\vartheta_f$	유속 자체를 음함수로
full-Newton	셋 모두	완전 전개

논문의 결과가 갈리는 지점이 여기다. 시간항에 Newton 선형화를 쓰는 것만으로 음파 전파 문제에서 1.4~1.5배 빨라진다. 이송항의 선택은 저마하에서는 거의 티가 나지 않다가, 마하 3 전방계단 문제에서 $Co = 0.9$ 로 시간스텝을 키우면 결정적으로 변한다. 그 조건에서 단일 루프 해법으로 수렴한 것은 $\rho$ -Newton 계열뿐이었다. $\vartheta$ -Newton을 더한 full-Newton은 수렴률의 음수 구간을 없애지만 실행 시간 이득은 작았다.

온도를 언제 갱신하는가 — 단일 루프와 이중 루프#

선형화만이 변수는 아니다. 비선형 반복의 구조도 갈린다.

단일 루프는 단순하다. 선형계를 풀고, 엔탈피에서 온도를 갱신하고, $p^{(n+1)}$ 과 $T^{(n+1)}$ 로 밀도를 갱신하고, 이송 속도를 갱신한 뒤 잔차를 본다.

\frac{\lVert A^{(n+1)}\varphi^{(n+1)} - b^{(n+1)}\rVert}{\lVert b^{(n+1)}\rVert\,\Theta} < \eta

$\Theta = \sqrt{N_r}$ 는 잔차 벡터 크기에 따른 정규화 인자다.

이중 루프는 Xiao 등의 처방을 따른다. 안쪽 루프에서는 밀도 갱신에 쓰는 온도를 상수로 고정한다. 즉 밀도를 압력만의 함수로 본다. 유동이 등온이라는 뜻이 아니다. 밀도를 계산할 때 쓰는 온도만 얼려두는 것이다. 안쪽이 수렴하면 바깥 루프에서 갱신된 온도로 밀도를 다시 계산한다. 이 구조가 언더릴랙세이션을 대신한다.

논문의 결론은 둘을 대립시키지 않는다. 모든 시간항과 이송항에 Newton 선형화를 일관되게 적용하면 언더릴랙세이션을 어떤 형태로도 쓰지 않을 수 있고, 그러면 단일 루프가 이중 루프보다 빨라진다. 마하 3 전방계단에서 단일 루프 $\rho$ -Newton은 10,616초, 같은 조건 이중 루프는 12,117초였다. 원뿔 초음속 유동에서도 순서가 같았다.

Python으로 좁혀보는 수렴 경계#

여기서부터는 논문의 코드가 아니라, 그 주장을 가장 작게 줄인 장난감 모형이다. 셀 하나, 등온 이상기체, 면 하나의 질량유속만 남긴다.

\rho(p) = \frac{p}{a_T^2}, \qquad \vartheta(p) = u_0 - \hat d\,(p - p_0)

$P = p/p_0$ 로 무차원화하고 목표 질량유속을 $\rho_0 u_0$ 로 두면 풀어야 할 비선형 방정식은 하나로 줄어든다.

m(P) = P\big[1 - D\,(P-1)\big] = 1, \qquad D = \frac{\hat d\,p_0}{u_0} = \frac{Co}{M}

여기에 두 선형화를 적용하면 서로 다른 반복 사상(map)이 나온다. 고정계수는

P^{(n+1)} = g_{\mathrm{fix}}(P^{(n)}) = 1 + \frac{1 - 1/P^{(n)}}{D}, \qquad \big|g_{\mathrm{fix}}'(1)\big| = \frac{1}{D} = \frac{M}{Co}

Newton 쪽은 정확히 $m(P)-1=0$ 에 대한 Newton–Raphson 반복과 일치한다. 대입하면 식은 그대로 떨어진다.

수축사상 조건 $|g'|<1$ 이 곧 $M < Co$ 다. 앞 절의 시소가 뒤집히는 지점과 같은 숫자다. 한 가지 더 있다. $m(P)=1$ 은 2차식이라 근이 둘이다. $P=1$ 과 $P=1/D=M/Co$ . 물리 근이 밀어내기 시작하면 반복은 발산하는 대신 다른 가지로 끌려가는 경우가 더 많다.

import math
 
def face_mass_flux(P, D):
    """비차원 면 질량유속 m/(rho0 u0).  P = p/p0,  D = Co/M."""
    return P * (1.0 - D * (P - 1.0))
 
def iterate_lagged(P0, D, nmax=40, eta=1e-10):
    """고정계수: rho^(n) theta^(n+1) = mdot*  — 밀도를 계수로 붙들어 둔다."""
    P, hist = P0, []
    for _ in range(nmax):
        r = abs(face_mass_flux(P, D) - 1.0)
        hist.append(r)
        if r < eta:
            return P, hist
        P = 1.0 + (1.0 - 1.0 / P) / D
        if not math.isfinite(P) or P <= 0.02 or P > 8.0:
            return float('nan'), hist
    return P, hist
 
def iterate_newton(P0, D, nmax=40, eta=1e-10):
    """Newton: rho^(n)theta^(n+1) + rho^(n+1)theta^(n) - rho^(n)theta^(n) = mdot*."""
    P, hist = P0, []
    for _ in range(nmax):
        r = abs(face_mass_flux(P, D) - 1.0)
        hist.append(r)
        if r < eta:
            return P, hist
        dm = (1.0 + D) - 2.0 * D * P           # d(rho theta)/dP
        P = P - (face_mass_flux(P, D) - 1.0) / dm
    return P, hist
 
def verdict(P, hist):
    if math.isnan(P):
        return "발산"
    if abs(P - 1.0) < 1e-4:
        return f"물리 근 ({len(hist)}회)"
    return f"다른 가지 P={P:.3f} ({len(hist)}회)"
 
cases = [("acoustic wave",  0.003, 0.10),
         ("Sod shock tube", 0.900, 0.40),
         ("forward step",   3.000, 0.90),
         ("forward step*",  3.000, 0.30)]
 
print(f"{'case':<16}{'M':>7}{'Co':>6}{'M/Co':>7}   {'고정계수':<22}{'Newton'}")
for name, M, Co in cases:
    D = Co / M
    Pl, hl = iterate_lagged(1.18, D)
    Pn, hn = iterate_newton(1.18, D)
    print(f"{name:<16}{M:>7.3f}{Co:>6.2f}{M / Co:>7.2f}   "
          f"{verdict(Pl, hl):<22}{verdict(Pn, hn)}")
 
print("\nSod shock tube — 반복별 잔차 ||r||")
_, hl = iterate_lagged(1.18, 0.40 / 0.90)
_, hn = iterate_newton(1.18, 0.40 / 0.90)
for n in range(5):
    print(f"  n={n}   고정계수 {hl[n]:.3e}   Newton {hn[n]:.3e}")

case                  M    Co   M/Co   고정계수                  Newton
acoustic wave     0.003  0.10   0.03   물리 근 (9회)             물리 근 (5회)
Sod shock tube    0.900  0.40   2.25   다른 가지 P=2.250 (32회)   물리 근 (5회)
forward step      3.000  0.90   3.33   다른 가지 P=3.333 (23회)   물리 근 (5회)
forward step*     3.000  0.30  10.00   발산                    물리 근 (4회)
 
Sod shock tube — 반복별 잔차 ||r||
  n=0   고정계수 8.560e-02   Newton 8.560e-02
  n=1   고정계수 1.383e-01   Newton 2.081e-02
  n=2   고정계수 1.725e-01   Newton 5.570e-04
  n=3   고정계수 1.565e-01   Newton 4.454e-07
  n=4   고정계수 1.061e-01   Newton 2.855e-13

Newton은 네 경우 모두 4~5회다. 잔차가 매 반복 제곱으로 줄어드는 2차 수렴이 그대로 보인다. 고정계수는 $M/Co$ 가 1을 넘는 순간 잔차가 세 번 연속 커진 뒤 다른 가지로 흘러간다. 위 네 칸을 마하 슬라이더와 fixed-coefficient/Newton 전환으로 그대로 따라가 보자.

linearisation

Mach M0.900Courant Co0.400initial guess P⁽⁰⁾1.180

paper test-cases

fixed-coefficient 상태에서 마하 슬라이더를 Courant 값 아래로 내리면 계단이 초록 점( $P=1$ )으로 조여든다. 반대로 올리면 같은 계단이 초록 점을 밀어내고 빨간 점( $P=M/Co$ )으로 걸어간다. Newton으로 바꾸면 어느 조합에서도 오른쪽 잔차 그래프가 서너 번 만에 $\eta$ 선 아래로 떨어진다.

장난감이라는 점은 분명히 해둔다. 실제 솔버에는 에너지 방정식, 다차원 이송, 압력 감쇠항의 시간 의존성이 더 얹힌다. $M = Co$ 는 정확한 경계가 아니라 눈금이다. 다만 논문이 보고한 순서 — 음파 문제는 어느 선형화로도 돌고, 전방계단은 $\rho$ -Newton 없이는 안 돈다 — 와 이 눈금이 같은 방향을 가리킨다.

시작값은 $P_0=1.18$ 이다. Sod의 $D=0.40/0.90=0.444$ 이면 $m(1.18)=1.18[1-0.444\times 0.18]=1.0856$ 이라 잔차 $|m-1|=8.56\times 10^{-2}$ 가 n=0에 그대로 찍힌다. 고정계수는 다음 두 번에 $1.383\times 10^{-1}$ , $1.725\times 10^{-1}$ 로 커진 뒤 $1.565\times 10^{-1}$ 에서 꺾인다. 줄어도 도착점은 $P=1$ 이 아니다. $P=1/D=M/Co=2.250$ 이다. 전방계단 $Co=0.9$ 도 같다. $3/0.9=3.333$ 이 23회 뒤의 다른 가지다. $Co=0.3$ 이면 $M/Co=10$ 이라 반복이 $P>8$ 에서 끊기고 발산으로 적힌다.

Newton 잔차는 $8.560\times 10^{-2}$ 에서 $2.081\times 10^{-2}$ , $5.570\times 10^{-4}$ , $4.454\times 10^{-7}$ , $2.855\times 10^{-13}$ 으로 떨어진다. $(2.081\times 10^{-2})^{2}\approx 4.3\times 10^{-4}$ 라서 다음 자리 $5.57\times 10^{-4}$ 와 자릿수가 같다. 음파만 고정계수도 물리 근이다. $M/Co=0.03$ 이라 $|g'(1)|$ 이 1보다 훨씬 작고, 9회 대 Newton 5회다. 논문의 마하 3 전방계단은 단일 루프 $\rho$ -Newton 10,616초, 이중 루프 12,117초로 약 1.14배다. 손잡이를 떼서 버는 시간은 그 정도고, 선형화를 잘못 고르면 그 전에 다른 근으로 간다.

M/Co가 1을 넘으면 고정계수는 다른 근에 앉는다#

Sod의 $M/Co=2.25$ , 전방계단의 3.33과 10은 모두 1보다 크다. 고정계수는 그 세 칸에서 물리 근 $P=1$ 을 내놓지 않았다. 32회와 23회 만에 $P=M/Co$ 로 가거나, 10에서는 발산했다. Newton은 같은 $P_0=1.18$ 에서 4회 또는 5회에 $P=1$ 이다. 언더릴랙세이션을 0.3으로 내리는 일은 이 분기를 고르지 않는다.

$M/Co=0.03$ 인 음파만 고정계수도 9회에 물리 근이다. 손잡이를 떼도 되는 칸과 떼면 다른 가지로 가는 칸이 같은 표에 있다. 시간스텝을 1001배 키운 반음해 압력 스킴이 $L_{\infty}$ 에서 흔들린 기록은, 반복 횟수가 줄었다고 그 스텝의 해가 안전한 것은 아니라는 쪽을 따로 보여 준다.