이완 계수 하나 바꿨더니 수렴해가 달라졌다 — Rhie–Chow 면속도가 언더릴랙세이션을 기억하는 자리

표준 Rhie–Chow는 수렴 후에도 압력 감쇠항 앞에 α_u를 남긴다. 면속도의 이전 값을 셀 평균이 아니라 면에 저장된 값으로 쓰면 수렴해가 α와 무관해진다.

이완 계수만 바꿨는데 수렴해가 달라졌다#

정상상태 SIMPLE 계산이 잔차 $10^{-8}$ 까지 내려갔다. 운동량 이완 계수 $\alpha_u$ 를 0.7에서 0.3으로 바꿔 다시 돌렸다. 반복은 세 배쯤 늘었지만 역시 수렴했다. 그런데 두 결과의 속도장이 다르다. 격자도, 경계조건도, 이산화도 같은데 말이다.

이완 계수는 반복 경로만 바꿔야 한다. 수렴점이 그 손잡이에 따라 움직이면 결과를 믿을 수 없다. 이 글은 그 원인이 Rhie–Chow 면속도 공식 한 줄에 있다는 것을 보인다. 1D 정수압 기둥 하나로 차이를 숫자로 재고, Majumdar(1988)의 보정이 왜 그 숫자를 0으로 만드는지 확인한다.

배경으로 SIMPLE과 Rhie–Chow의 기본 구조를 알고 있으면 충분하다. 동위치 격자(collocated grid, 속도와 압력을 같은 셀 중심에 두는 배치)를 전제로 한다.

직접 돌려보는 두 층 기둥#

α_u 슬라이더와 standard/Majumdar 토글이, 같은 기둥에서 셀 속도 막대가 노란 눈금에서 얼마나 밀리는지를 바꾼다.

alpha_u = 0.30 density ratio = 10 iters/frame = 2

iterating | iter = 0 | residual = 1.0e+0 | peak |u_P| = 0.0000 | cells |u|>0.01 = 0 | max |p-p(a=1)| = 0.00e+0

standard에서 $\alpha_u$ 를 0.1까지 내리면 막대 부호가 칸마다 바뀌며 넓게 번진다. Majumdar로 바꾸면 어떤 $\alpha_u$ 에서도 막대가 눈금 위로 돌아온다.

SIMPLE 한 바퀴에서 α가 들어가는 자리#

셀 P의 이산 운동량 방정식은 다음 꼴이다.

a_P u_P = \sum_{nb} a_{nb} u_{nb} + b_P - V_P (\nabla p)_P

$a_P$ 는 대각 계수, $a_{nb}$ 는 이웃 계수, $b_P$ 는 압력 이외의 소스, $V_P$ 는 셀 부피다.

SIMPLE은 이 식을 그대로 풀지 않는다. 대각을 $a_P/\alpha_u$ 로 키우고, 그만큼을 이전 반복값 $u_P^{old}$ 로 오른쪽에 보탠다.

\frac{a_P}{\alpha_u} u_P = \sum_{nb} a_{nb} u_{nb} + b_P - V_P (\nabla p)_P + \frac{1-\alpha_u}{\alpha_u} a_P u_P^{old}

수렴하면 $u_P = u_P^{old}$ 이므로 추가한 두 항은 상쇄된다. 셀 중심 운동량 방정식만 놓고 보면 이완은 수렴점에 흔적을 남기지 않는다. 문제는 면에서 생긴다.

면속도 공식에 $u^{old}$ 가 섞이는 경로#

이완된 식을 셀 속도에 대해 다시 쓰면 다음과 같다.

u_P = \hat{u}_P - d_P (\nabla p)_P,\qquad d_P = \frac{\alpha_u V_P}{a_P}

\hat{u}_P = \frac{\alpha_u}{a_P}\Big(\sum_{nb} a_{nb} u_{nb} + b_P\Big) + (1-\alpha_u)\,u_P^{old}

$\hat{u}_P$ 는 압력 기울기를 뺀 나머지 속도, $d_P$ 는 압력 기울기에 대한 속도 민감도다.

Rhie–Chow 보간은 면 f에서 $\hat{u}$ 와 $d$ 는 보간하고, 압력 기울기만 면의 두 셀 값으로 다시 계산한다.

u_f = \overline{\hat{u}}_f - d_f (\nabla p)_f

$\overline{\,\cdot\,}_f$ 는 두 셀의 선형 보간, $(\nabla p)_f = (p_N - p_P)/\Delta x$ 는 콤팩트한 면 기울기다.

문제는 $\overline{\hat{u}}_f$ 안에 $(1-\alpha_u)\,\overline{u^{old}}_f$ 가 들어 있다는 점이다. 면속도의 "이전 값"을 셀 속도의 평균으로 대신 쓴 것이다. 셀 속도 평균은 면속도와 같지 않다. 그 차이가 바로 Rhie–Chow 감쇠항이다.

수렴하면 α가 감쇠항 앞에 남는다#

수렴 상태에서 $u^{old} = u$ 를 넣으면 표준 공식은 다음이 된다.

u_f = \overline{u}_f + \alpha_u \,\overline{\Big(\tfrac{V}{a}\Big)}_f \Big[\,\overline{(\nabla p)}_f - (\nabla p)_f\,\Big]

대괄호는 셀 기울기 평균과 콤팩트 면 기울기의 차이다. 매끄러운 압력에서는 $O(\Delta x^2)$ 이고, 체커보드(checkerboard, 격자 간격으로 부호가 바뀌는 압력 모드)에서는 가장 크다. 이 항이 체커보드를 연속 방정식에 보이게 만드는 장치다.

이 장치의 세기가 $\alpha_u$ 배로 줄어든다. $\alpha_u = 0.1$ 이면 감쇠는 10%만 남는다. 수렴해가 $\alpha_u$ 에 의존하고, $\alpha_u$ 를 줄일수록 simpleFoam에서 압력 잔차가 안 죽는 격자 무늬가 슬그머니 되돌아온다.

Majumdar의 처방은 간단하다. 면속도의 이전 값을 면에 따로 저장해 그대로 쓴다.

u_f = \overline{\hat{u}^{*}}_f - d_f (\nabla p)_f + (1-\alpha_u)\, u_f^{old}

$\hat{u}^{*}$ 는 $\hat{u}$ 에서 $(1-\alpha_u)u^{old}$ 항을 뺀 값이다. 수렴하면 $u_f = u_f^{old}$ 이므로 $(1-\alpha_u)$ 가 양변에서 정리되고, 감쇠항 앞의 계수는 1로 돌아간다.

코드: 같은 기둥, 다섯 개의 α#

시험 문제는 일부러 정답이 자명한 것을 골랐다. 양끝이 벽인 1D 기둥에 아래 절반은 $\rho = 10$ , 위 절반은 $\rho = 1$ 을 채우고 중력을 건다. 정답은 정지 상태 $u = 0$ 이다. 1D에서는 연속 방정식이 모든 면에서 $u_f = 0$ 을 강제하므로, 압력 보정은 적분 한 번으로 끝난다. 남는 셀 속도 $u_P$ 가 곧 보간 오차다.

# 1D 정수압 기둥: 아래 절반 rho=10, 위 절반 rho=1, 양끝 벽. 정답은 u = 0.
def density_column(n, rho_lo=1.0, rho_hi=10.0):
    return [rho_hi if (i + 0.5) / n < 0.5 else rho_lo for i in range(n)]
 
def cell_grad(p, h):
    # 면 압력 = 이웃 평균, 벽 면은 선형 외삽
    n = len(p)
    pf = [1.5 * p[0] - 0.5 * p[1]]
    pf += [0.5 * (p[i] + p[i + 1]) for i in range(n - 1)]
    pf += [1.5 * p[-1] - 0.5 * p[-2]]
    return [(pf[i + 1] - pf[i]) / h for i in range(n)]
 
def simple_hydrostatic(n, alpha, majumdar, a=1.0, g=1.0, tol=1e-13):
    h = 1.0 / n
    f = [-r * g for r in density_column(n)]
    p, u, uf = [0.0] * n, [0.0] * n, [0.0] * (n + 1)
    d = alpha / a                      # d = alpha*V/a_P (V = 1 단위)
    for it in range(1, 5000):
        gp = cell_grad(p, h)
        # 이완된 운동량 예측자
        us = [alpha * (f[i] - gp[i]) / a + (1 - alpha) * u[i] for i in range(n)]
        # Rhie-Chow: 압력항을 뺀 u_hat을 면으로 보간
        uh = [us[i] + d * gp[i] - majumdar * (1 - alpha) * u[i] for i in range(n)]
        ufs = [0.0] * (n + 1)
        for i in range(1, n):
            ufs[i] = 0.5 * (uh[i - 1] + uh[i]) - d * (p[i] - p[i - 1]) / h
            ufs[i] += majumdar * (1 - alpha) * uf[i]
        # 압력 보정: 1D에서는 연속식이 면속도 0을 바로 요구한다
        pc = [0.0] * n
        for i in range(1, n):
            pc[i] = pc[i - 1] + h * ufs[i] / d
        p = [p[i] + pc[i] for i in range(n)]
        gc = cell_grad(pc, h)
        un = [us[i] - d * gc[i] for i in range(n)]
        uf = [0.0] * (n + 1)          # 보정 후 면속도는 정확히 0
        res = max(abs(un[i] - u[i]) for i in range(n))
        u = un
        if res < tol and it > 5:
            return u, p, it
    return u, p, it
 
n = 80
_, p_ref, _ = simple_hydrostatic(n, 1.0, False)
print("alpha | iters | peak|u|  cells>0.01 | max|p-p(1)|   <- standard / majumdar")
for alpha in (1.0, 0.7, 0.5, 0.3, 0.1):
    row = []
    for mj in (0, 1):
        u, p, it = simple_hydrostatic(n, alpha, mj)
        wide = sum(abs(v) > 0.01 for v in u)
        dp = max(abs(p[i] - p_ref[i]) for i in range(n))
        row.append(f"{it:4d} {max(map(abs, u)):.4f} {wide:3d} {dp:.2e}")
    print(f" {alpha:.1f}  | " + "  /  ".join(row))
 
u, _, _ = simple_hydrostatic(n, 0.1, 0)
print("standard, alpha=0.1, cells 34..45:", " ".join(f"{v:+.3f}" for v in u[34:46]))

alpha | iters | peak|u|  cells>0.01 | max|p-p(1)|   <- standard / majumdar
 1.0  |    6 2.2500   2 0.00e+00  /     6 2.2500   2 0.00e+00
 0.7  |   26 2.0499   6 4.59e-03  /    27 2.2500   2 8.88e-16
 0.5  |   42 1.8640   6 8.24e-03  /    45 2.2500   2 0.00e+00
 0.3  |   79 1.5925  10 1.27e-02  /    84 2.2500   2 8.88e-16
 0.1  |  247 1.0811  16 1.92e-02  /   271 2.2500   2 0.00e+00
standard, alpha=0.1, cells 34..45: +0.041 -0.079 +0.152 -0.292 +0.562 -1.081 +1.081 -0.562 +0.292 -0.152 +0.079 -0.041

1.0811과 16칸이 말하는 것#

$\alpha_u = 1$ 에서는 두 공식이 같다. 계면 양쪽 두 셀에 $\mp 2.25$ 가 남는다. 이것은 이완과 무관한 다른 오차다. 체적력 $\rho g$ 의 계단을 Rhie–Chow 면이 보지 못해서 생긴다. 밀도비 $r$ 에서 크기는 $(r-1)g/4$ 이고, 위젯의 밀도비 슬라이더로 확인할 수 있다. Denner(2013) 같은 balanced-force 처리가 이 항을 면에서 압력과 같은 방식으로 다뤄 지운다.

표준 공식은 $\alpha_u$ 를 낮출수록 최대 속도가 2.25에서 1.0811로 줄어든다. 숫자만 보면 좋아진 것 같다. 그러나 $|u| > 0.01$ 인 셀은 2개에서 16개로 늘었다. 마지막 줄을 보면 부호가 칸마다 바뀌고, 진폭이 한 칸에 약 1.9배씩 줄어든다. 감쇠가 약해진 만큼 오차가 국소에 머물지 못하고 진동 모드로 퍼진 것이다.

압력도 움직였다. $\alpha_u = 0.1$ 에서 $\alpha_u = 1$ 해와의 차이는 $1.92\times10^{-2}$ 다. Majumdar 쪽은 다섯 경우 모두 $10^{-15}$ 이하로, 기계 정밀도에서 같은 해다. 반복 수는 Majumdar가 약간 더 많지만 수렴 경로의 차이일 뿐이다.

반복 수를 공식끼리 나란히 놓으면, 이완이 경로만 늘리는 일과 해를 바꾸는 일이 갈린다. $\alpha_u = 1$ 이면 둘 다 6번이다. 0.7이면 26과 27, 0.5면 42와 45, 0.3이면 79와 84, 0.1이면 247과 271이다. 0.7에서 0.3으로 내리면 표준 쪽은 $79/26 \approx 3.0$ 배다. 앞에서 말한 세 배가 이 칸이다. Majumdar는 한 박자 늦을 뿐, 비율은 비슷하다.

해가 갈리는 열은 peak $|u|$ 다. 표준은 2.2500, 2.0499, 1.8640, 1.5925, 1.0811로 내려간다. $\alpha_u = 0.1$ 이면 감쇠항 앞의 계수가 0.1이라 체커보드를 누르는 힘이 10%만 남는다. 그런데 최대 속도는 2.25의 10%인 0.225가 아니라 1.0811이다. 오차가 준 것이 아니다. 계면의 $\mp 2.25$ 가 이웃으로 번졌다. $|u| > 0.01$ 인 셀이 2, 6, 6, 10, 16으로 늘어난 것이 그 번짐이다.

34번부터 부호가 바뀌는 열의 진폭 비는 거의 일정하다. $1.081/0.562 \approx 1.92$ , $0.562/0.292 \approx 1.92$ , $0.292/0.152 \approx 1.92$ 다. 한 칸에 약 1.9배. 압력 차도 $4.59\times10^{-3}$ , $8.24\times10^{-3}$ , $1.27\times10^{-2}$ , $1.92\times10^{-2}$ 로 $\alpha_u$ 가 작을수록 커진다. Majumdar 열은 다섯 $\alpha$ 모두 peak $|u| = 2.2500$ , 넓은 셀은 2개, 압력 차는 0 또는 $8.88\times10^{-16}$ 이다. 2.25는 $(r-1)g/4$ 에 $r = 10$ , $g = 1$ 을 넣은 값이다. $\alpha_u$ 는 반복 횟수에만 남는다.

Majumdar 보정이 저장하는 것#

구현 비용은 면마다 실수 하나다. 압력 보정이 끝난 뒤 보정된 면 플럭스를 $u_f^{old}$ 로 저장하면 된다. 대부분의 코드는 이미 질량 플럭스를 면에 보관하므로 사실상 공짜다.

주의할 점이 하나 있다. 저장하는 값은 보간 직후의 $u_f^{*}$ 가 아니라 압력 보정까지 마친 면속도다. 보정 전 값을 저장하면 연속 방정식을 만족하지 않는 플럭스가 다음 반복의 기억이 된다.

시간항도 같은 함정을 판다#

비정상 계산에서는 시간항이 같은 역할을 한다. 음해 시간 이산화의 $\rho V u_P^n/\Delta t$ 가 $\hat{u}_P$ 에 섞여 셀 평균으로 보간되면, 정상상태에 도달한 해가 $\Delta t$ 에 따라 달라진다. Choi(1999) 등이 이 의존성을 지적했다. OpenFOAM의 ddtCorr는 이전 시간 단계의 면 플럭스를 직접 써서 이 차이를 보정하는 항이다.

원리는 이완과 똑같다. "이전 값"이 필요한 모든 자리에서 면에는 면의 값을 쓴다. 여러 시간 단계를 쓰면 저장할 면 플럭스도 그만큼 늘어난다. 밀도 계단이 남기는 $\mp 2.25$ 처럼 정지해가 흔들리는 자리는 극좌표에서 well-balanced가 깨질 때와 같은 종류의 균형 오차다.

이완 계수를 수렴 속도 손잡이로만 남기려면#

검사는 간단하다. 같은 케이스를 $\alpha_u$ 두 값으로 수렴시키고 결과를 빼 본다. 차이가 반복 허용오차 수준이면 면속도 공식은 깨끗하다. 차이가 격자 오차와 비슷한 크기라면 면 보간이 $u^{old}$ 를 셀 평균으로 섞고 있을 가능성이 높다.

비정상 해석이라면 같은 검사를 $\Delta t$ 로 반복한다. 정상상태로 수렴한 해가 $\Delta t$ 에 따라 달라지면 시간항 쪽 보정이 빠진 것이다. 두 검사 모두 코드 한 줄을 의심하게 하는 데 충분하다.