임의 방향에서 풀어낸 오일러 고유 구조 — 세 파동, 두 중복, 한 특이점

출장지 호텔에서 노트북을 열었다. 사흘 전 격자를 임의 방향으로 회전시킨 뒤로 솔버가 자꾸 음의 밀도를 뱉었다. $n_x=1, n_y=0$ 인 정렬 케이스에선 멀쩡했다. $n_x=0.01$ 만 되어도 폭주가 시작됐다. 로그를 거꾸로 더듬으니 한 줄의 부호가 갈렸다. 그 줄은 좌측 고유벡터 행렬의 마지막 행이었다. 사흘 동안 이 한 줄이 무엇을 숨기고 있는지를 알아내야 했다.

이 글은 그 사흘의 정리다. 임의 단위 법선 $\hat n = (n_x, n_y, n_z)$ 위에서 3D Euler 방정식의 고유값과 좌·우 고유벡터를 한 번에 적고, 좌측 고유벡터가 왜 한 컨벤션으로는 항상 무너지는지를 보인다. 끝에서 같은 식으로 Roe-type 플럭스를 짠다.

보존형 오일러를 임의 방향으로 본다#

3차원 비점성 오일러 방정식은 보존형 적분식이다.

여기서 $\mathbf{Q} = (\rho,\rho u,\rho v,\rho w,\rho e_o)^T$ 는 보존 변수, $\mathbf{F}\cdot \hat n$ 은 법선 방향 플럭스다. 셀 면 위에서 우리가 다루는 양은 한 면에 묶여 있는 법선 속도 $v_n = u n_x + v n_y + w n_z$ 다. 음속 $a$ 는 $a^2 = \gamma R T = \gamma p/\rho$. 단위 법선이라 $n_x^2 + n_y^2 + n_z^2 = 1$.

이 좌표가 좋은 이유는 단순하다. 면에 수직인 한 변수 $v_n$ 으로 모든 파동이 정리되기 때문이다.

야코비안 A(n̂) — 5×5 안에 들어 있는 다섯 파동#

$\mathbf{F}\cdot \hat n$ 을 $\mathbf{Q}$ 로 미분하면 변환 행렬 $A(\hat n)= \partial(\mathbf F\cdot\hat n)/\partial\mathbf Q$ 를 얻는다. 이 행렬은 전역 $(x,y,z)$ 나 국소 $(\xi,\eta,\zeta)$ 로 따로 적을 필요가 없다. 한 번에 $\hat n$ 만 들고 적으면 된다. 이것이 1983년 Yee–Kutler 이후로 굳어진 표준이다.

특성방정식 $\det(A - \lambda I) = 0$ 의 근은 다섯 개다.

여기서 $v_n$ 은 세 번 나타나는데, 두 번은 중복근이다. 이 중복이 좌측 고유벡터를 흔드는 첫 번째 단서다.

위 그림에서 슬라이더로 $M = v_n/a$ 를 흔들어 보면, 세 광선이 부채꼴로 벌어진다. $|M|<1$ 일 때는 원점을 사이에 두고 음향파가 양쪽으로 갈리고, $|M|>1$ 이 되는 순간 모든 광선이 같은 부호의 $x$ 쪽으로 휘어진다. 입구·출구 경계조건이 마하수에 따라 갯수를 바꾸는 이유가 이 한 장에 들어 있다.

세 가족, 두 중복 — 음향파와 엔트로피파#

다섯 고유값은 세 가족이다. 흐름 속도 $v_n$ 자체는 엔트로피파(entropy wave)이고, 이 파동은 엔트로피와 접선 속도라는 두 가지 별개의 신호를 같은 속도로 운반한다. 그래서 중복이 둘이다.

$v_n \pm a$ 는 음향파(acoustic wave) 두 갈래로, 이쪽은 압력 신호와 법선 속도 신호를 매단다. Riemann 문제의 부채꼴 그림에서 가운데 직선이 접촉(contact), 양옆의 부채/충격이 음향파다.

이 구조를 손으로 외워두면 코드 디버깅이 빨라진다. 밀도만 점프가 있고 $\Delta p = \Delta u = 0$ 이라면 엔트로피파 강도만 살아남고 양쪽 음향파 강도는 0 이어야 한다. 이걸 한 번에 확인할 그림이 아래다.

"no jump" 프리셋에서 시작해 $\rho_R$ 만 조금씩 흔들어보자. 가운데 회색 막대만 자라고 양쪽 빨강·파랑은 0 에 붙어 있다. 이번엔 "Sod" 로 바꾸자. 우측 음향파 막대가 크게 자란다 — 그 자리가 충격이다. "123 problem" 은 좌우 대칭 팽창파라서 두 음향 막대가 같은 부호로 크게 뜨고 접촉은 거의 0 이다.

우측 고유벡터 — 한 번에 적기#

$A R_i = \lambda_i R_i$ 의 해는 5×5 행렬로 묶인다. 한 컨벤션(R-1)으로 적으면

각 기호 해석은 다음 한 줄이다. $e_k = \tfrac12(u^2+v^2+w^2)$ 는 단위 질량당 운동에너지, $h_o = h + e_k$ 는 단위 질량당 정체엔탈피, 첫 세 열은 각각 $\lambda_1, \lambda_2, \lambda_3$ 에 대응한다. 마지막 두 열은 중복근 $\lambda_4=\lambda_5=v_n$ 에 대응하는 두 접선 방향 벡터이며, 컨벤션마다 자리가 다르다.

세 번째와 다섯 번째 열이 "어떤 접선 방향을 골랐느냐"의 흔적이다. 다른 두 컨벤션(R-2, R-3)은 같은 부분 공간을 다른 두 기저로 펼친 것일 뿐이다.

좌측 고유벡터의 함정 — n_x → 0 이 부수는 줄#

$L = R^{-1}$ 을 구해보면 마지막 두 행에 $1/n_x$ 가 등장한다. R-1 컨벤션에선 항상 $n_x$ 로 나누는 줄이 있다. $\hat n = (0, 1, 0)$ 처럼 $n_x = 0$ 인 면이 한 번이라도 등장하면 그 면의 좌측 고유벡터는 발산한다. 이게 출장지에서 사흘 잡아먹은 버그였다.

해결은 단순하다. R-2 는 $n_y$ 로, R-3 는 $n_z$ 로 나누게 만들어두면 된다. 각 면마다 $|n_x|, |n_y|, |n_z|$ 중 가장 큰 성분을 분모로 잡는 컨벤션을 고르면 된다. 코드 한 줄로 표현하면 다음과 같다.

def pick_convention(n):
    # n: ndarray of shape (3,), unit normal
    k = int(np.argmax(np.abs(n)))  # 0, 1, or 2
    return k  # use R-1, R-2, or R-3 accordingly

세 컨벤션 중 하나는 그 면에서 반드시 안전하다는 게 핵심 결과다. 비정렬 격자 코드는 이 한 줄을 면별로 호출하지 않으면 살아남지 못한다.

2D로 좁히면 사라지는 특이점#

$w=n_z=0$ 으로 줄이면 5×5 가 4×4 가 된다. 다섯 고유값은 $\{v_n - a, v_n, v_n + a, v_n\}$ 네 개로 줄고, 좌측 고유벡터 행렬의 마지막 행은 모두 두 항이 만나 깔끔히 상쇄된다. 예를 들어

처럼, 형식상 분모에 $n_x$ 가 있어도 분자가 $n_x$ 로 나누어진다. 그래서 2D 코드만 가지고 살던 사람은 이 함정을 한 번도 못 만나고 평생을 보낼 수 있다. 3D 로 옮기는 순간 새로 나타나는 한 줄이다.

Python — 임의 방향 Roe-type 플럭스 한 함수로#

세 가족 고유값과 1D Roe-평균 파동 강도를 한 곳에 모아, $\hat n$ 방향 Roe-type 플럭스를 짠다. 토이 문제는 비스듬한 충격관: $\hat n = (\cos 30°, \sin 30°, 0)$ 면에 양쪽 상태를 던진다.

import numpy as np
 
GAMMA = 1.4
 
def primitive_to_state(rho, vel, p):
    """rho: scalar, vel: (3,), p: scalar -> conservative Q (5,)"""
    rE = p / (GAMMA - 1) + 0.5 * rho * np.dot(vel, vel)
    return np.array([rho, rho * vel[0], rho * vel[1], rho * vel[2], rE])
 
def flux_normal(Q, n_hat):
    rho = Q[0]
    vel = Q[1:4] / rho
    rE = Q[4]
    p = (GAMMA - 1) * (rE - 0.5 * rho * np.dot(vel, vel))
    vn = np.dot(vel, n_hat)
    H = (rE + p) / rho
    return np.array([
        rho * vn,
        rho * vel[0] * vn + p * n_hat[0],
        rho * vel[1] * vn + p * n_hat[1],
        rho * vel[2] * vn + p * n_hat[2],
        rho * H * vn,
    ])
 
def roe_normal_flux(QL, QR, n_hat):
    rL, rR = QL[0], QR[0]
    vL = QL[1:4] / rL
    vR = QR[1:4] / rR
    pL = (GAMMA - 1) * (QL[4] - 0.5 * rL * np.dot(vL, vL))
    pR = (GAMMA - 1) * (QR[4] - 0.5 * rR * np.dot(vR, vR))
    HL = (QL[4] + pL) / rL
    HR = (QR[4] + pR) / rR
 
    sL, sR = np.sqrt(rL), np.sqrt(rR)
    inv = 1.0 / (sL + sR)
    rho = sL * sR
    vel = (sL * vL + sR * vR) * inv
    H   = (sL * HL + sR * HR) * inv
    vn  = np.dot(vel, n_hat)
    a2  = (GAMMA - 1) * (H - 0.5 * np.dot(vel, vel))
    a   = np.sqrt(max(a2, 1e-12))
 
    eigs = np.array([vn - a, vn, vn + a])
    dvel = vR - vL
    dvn  = np.dot(dvel, n_hat)
    drho = rR - rL
    dp   = pR - pL
 
    alpha_contact = drho - dp / a2
    alpha_minus   = (dp - rho * a * dvn) / (2 * a2)
    alpha_plus    = (dp + rho * a * dvn) / (2 * a2)
 
    # average flux + |lambda| * alpha * R
    Favg = 0.5 * (flux_normal(QL, n_hat) + flux_normal(QR, n_hat))
 
    # build the three relevant right eigenvectors on the fly
    R_minus = np.array([1, vel[0] - a * n_hat[0], vel[1] - a * n_hat[1],
                        vel[2] - a * n_hat[2], H - a * vn])
    R_cont  = np.array([1, vel[0], vel[1], vel[2], 0.5 * np.dot(vel, vel)])
    R_plus  = np.array([1, vel[0] + a * n_hat[0], vel[1] + a * n_hat[1],
                        vel[2] + a * n_hat[2], H + a * vn])
 
    diss = (np.abs(eigs[0]) * alpha_minus   * R_minus
          + np.abs(eigs[1]) * alpha_contact * R_cont
          + np.abs(eigs[2]) * alpha_plus    * R_plus)
 
    # tangential velocity jump (carried by contact, lambda = vn)
    dvel_tan = dvel - dvn * n_hat
    R_tan = np.zeros(5); R_tan[1:4] = dvel_tan
    R_tan[4] = rho * np.dot(vel, dvel_tan)
    diss = diss + np.abs(vn) * R_tan
 
    return Favg - 0.5 * diss
 
if __name__ == '__main__':
    n_hat = np.array([np.cos(np.pi/6), np.sin(np.pi/6), 0.0])
    QL = primitive_to_state(1.0,   np.array([0.0, 0.0, 0.0]), 1.0)
    QR = primitive_to_state(0.125, np.array([0.0, 0.0, 0.0]), 0.1)
    F = roe_normal_flux(QL, QR, n_hat)
    print('Sod tilted by 30 deg, Roe flux =', F)

손으로 짠 형태라 빠르진 않다. 그러나 면 하나에 들어가는 식이 모두 노출되어 있어 디버깅에 좋다. 실제 운영용 코드는 pick_convention 으로 좌측 고유벡터 컨벤션을 면별 분기하고, 엔트로피 픽스(Harten 또는 LLF 부분 적용)를 마지막에 끼워야 한다.

np.argmax(np.abs(n_hat)) 한 줄을 잊은 채 R-1 한 컨벤션만 박아두면 — 사흘 잃는다. 출장지에서 잃은 건 호텔 조식 값보다 그 사흘이었다.

디버깅 일기를 닫으며#

• 임의 단위 법선 $\hat n$ 으로 적은 5×5 야코비안은 한 식으로 전역 $(x,y,z)$ 와 국소 $(\xi,\eta,\zeta)$ 를 통일한다.
• 고유값은 세 가족 — 음향파 $v_n \pm a$, 엔트로피파 $v_n$ 이중복. 중복이 좌측 고유벡터에 부분 공간을 만든다.
• 좌측 고유벡터 R-1 컨벤션은 $n_x=0$ 면에서 발산한다. 면마다 $\max |n_i|$ 인덱스로 컨벤션을 골라잡는 한 줄이 비정렬 코드를 살린다.