출구가 음파를 반사하지 않게 하는 법 — Navier–Stokes 특성 경계조건(NSCBC)

압축성 LES가 6시간 만에 발산했다. 로그에는 NaN도, 충격파도 없었다. 평균 압력이 천천히 0.3 bar 떠올랐고, 음파가 출구에서 도메인 안으로 튕겨 들어왔다. 범인은 출구의 zero-gradient extrapolation이었다. 오늘은 그 출구를 음파에게 "투명"하게 만드는 법 — Poinsot–Lele(1992)의 NSCBC를 단 하나의 1D 음파 펄스로 끝까지 따라간다.

6시간 만에 도메인이 발산했다#

압축성 솔버의 외부 경계는 충격파보다 더 자주 시뮬레이션을 죽인다. zero-gradient 출구는 운동량과 에너지를 "복사"해 ghost 셀에 붙이지만, 그 복사는 음파의 임피던스 부정합을 만든다. 도메인을 빠져나가야 할 음파가 임피던스 점프를 만나면, 그 일부는 반드시 안으로 되튄다.

특히 LES·DNS·연소 해석에서는 이 반사가 평균장에 누적된다. Reynolds 평균이 잡히기도 전에 출구가 "공명기" 역할을 시작하면, 압력이 천천히 표류(drift)한다. NSCBC는 이 문제를 "들어오는 특성파 진폭을 직접 지정한다"는 한 가지 원리로 풀어낸다.

출구가 보낸 다섯 신호#

1D Euler 방정식을 보존변수 $U = (\rho, \rho u, E)^T$ 에 대해 쓰면,

\frac{\partial U}{\partial t} + A(U)\frac{\partial U}{\partial x} = 0

여기서 $A = \partial F/\partial U$ 는 플럭스 야코비안이다. 이 야코비안의 우유사 분해는 다섯 개(3D 평면 경계의 경우)의 고유값 $\lambda_i \in \{u-c,\, u,\, u,\, u,\, u+c\}$ 을 준다. $u$ 는 유속, $c$ 는 음속이다.

각 고유값은 특성파(characteristic wave) 하나에 대응한다. $u-c$ 는 음향파(왼쪽 진행), $u$ 는 엔트로피·접선속도 정보(유체와 함께 이동), $u+c$ 는 음향파(오른쪽 진행)다.

Mach number M0.30

위 그림에서 Mach 수를 0에서 1.6까지 움직여 보자. 아음속 출구는 $u-c$ 하나만 안으로 들어온다(굵은 빨간선). 즉 BC가 필요한 변수는 정확히 한 개. 이걸 모르면 압력·속도·온도를 한꺼번에 강제하다가 over-specified가 되어 솔버가 죽는다.

LODI — 평면파만 통과시키는 단순화#

NSCBC의 두뇌는 Local One-Dimensional Inviscid(LODI) relations다. 경계에서 점성과 접선 미분이 한 시간 스텝 동안 무시 가능하다고 가정하면, 보존방정식은 특성 진폭 $\mathcal{L}_i$ 의 ODE가 된다.

\frac{\partial \rho}{\partial t} + \frac{1}{c^2}\left[\mathcal{L}_2 + \tfrac{1}{2}(\mathcal{L}_5 + \mathcal{L}_1)\right] = 0

\frac{\partial u}{\partial t} + \frac{1}{2\rho c}(\mathcal{L}_5 - \mathcal{L}_1) = 0

\frac{\partial p}{\partial t} + \frac{1}{2}(\mathcal{L}_5 + \mathcal{L}_1) = 0

여기서 $\mathcal{L}_1 = (u-c)\left(\partial p/\partial x - \rho c\, \partial u/\partial x\right)$ , $\mathcal{L}_5 = (u+c)\left(\partial p/\partial x + \rho c\, \partial u/\partial x\right)$ , $\mathcal{L}_2$ 는 엔트로피파이다.

핵심: 출구로 나가는 파 $\mathcal{L}_5,\,\mathcal{L}_2$ 는 내부 셀에서 upwind 차분으로 직접 계산할 수 있다. 들어오는 파 $\mathcal{L}_1$ 만 BC가 정해야 한다. 이 한 줄이 NSCBC 전부다.

비반사의 함정 — drifting mean pressure#

가장 단순한 비반사 BC는 들어오는 파를 0으로 둔다.

\mathcal{L}_1 = 0

음향파가 깨끗하게 빠져나간다. 그러나 LODI를 적분해 보면 $\partial p/\partial t = -\mathcal{L}_5/2$ 로, $\mathcal{L}_5$ 의 시간 평균이 작은 양수면 $p$ 가 무한정 표류한다. 코드는 안 죽지만, 6시간 뒤 평균 압력이 슬그머니 1 atm을 넘는다.

Rudy & Strikwerda(1980)와 Poinsot & Lele(1992)는 이 표류를 잡는 1차 완화 항을 제안했다.

\mathcal{L}_1 = K\,(p - p_\infty), \qquad K = \frac{\sigma\,(1 - M_{\max}^2)\,c}{L}

$p_\infty$ 는 타겟 압력, $L$ 은 도메인 특성 길이, $M_{\max}$ 는 경계에서 예상되는 최대 Mach. $\sigma$ 는 완화 강도. 추천값은 $\sigma \approx 0.25$ .

Poinsot–Lele 완화 — σ 한 변수로 잡는다#

$\sigma$ 가 가진 두 얼굴은 다음과 같다.

$\sigma$ 가 너무 작으면 ( $\sim 0$ ): 거의 완전 비반사. 평균 표류는 통제 불가.
$\sigma$ 가 너무 크면 ( $\sim 5$ ): 경계가 다시 stiff해진다. 음파가 다시 반사된다.

$\sigma$ 가 acoustic reflection coefficient $R$ 과 가지는 관계는 1D 선형 해석으로 $|R| \approx \sigma L/(c\tau)$ 정도로 근사된다 ( $\tau$ 는 파동 통과 시간). 따라서 LES에서는 도메인 길이와 평균 Mach가 정해진 뒤 $\sigma$ 를 0.1–0.5 사이에서 튜닝한다.

sigma (reflectivity)0.05

아래 시뮬레이션에서 직접 조작해 보자. 모드를 "Zero-gradient"에 두면 Gaussian 음파 펄스가 오른쪽 경계에서 거의 그대로 튕긴다 (반사 계수 ≈ 1). "NSCBC"로 바꾸면 펄스가 깔끔히 사라지지만, $\sigma$ 를 0.5 이상 올리면 잔향이 다시 보인다. 하단 그래프는 경계에서의 $p'$ 시간 이력이다.

Python — 음파 펄스가 출구를 만날 때#

리만 변수 두 개 ( $\mathcal{A}^+ = u + p/(\rho c)$ , $\mathcal{A}^- = u - p/(\rho c)$ )로 전이된 1D 선형화 Euler 모델이다. 토이 문제: 정지한 매질에 우향 가우시안 펄스 하나를 던지고 두 BC를 비교한다.

import numpy as np
 
def linearized_acoustic_pulse(bc='nscbc', sigma=0.05, N=240, n_steps=400):
    """1D linearized Euler with right-side acoustic outflow.
    Returns p'(x) at final time for the two BC modes."""
    # 격자 / 음속 / CFL
    c, cfl = 1.0, 0.5
    dx = 1.0 / N
 
    Ap = np.exp(-((np.linspace(0, 1, N) - 0.22) ** 2) / 0.003)  # 우향 펄스
    Am = np.zeros(N)                                             # 좌향 0
 
    for _ in range(n_steps):
        Ap_new = Ap.copy()
        Am_new = Am.copy()
 
        # Ap 는 +c 로 이동 → 좌측 upwind
        Ap_new[1:] = Ap[1:] - cfl * (Ap[1:] - Ap[:-1])
        Ap_new[0] = 0.0  # 유입 정지
 
        # Am 는 -c 로 이동 → 우측 upwind
        Am_new[:-1] = Am[:-1] - cfl * (Am[:-1] - Am[1:])
 
        # 우측 경계 — 핵심
        if bc == 'wall':
            # zero-gradient: 들어오는 Am 가 Ap 와 같다 → 완전 반사
            Am_new[-1] = Ap_new[-1]
        else:
            # NSCBC: 들어오는 Am 를 반사계수 σ 로 묶는다
            Am_new[-1] = sigma * Ap_new[-1]
 
        Ap, Am = Ap_new, Am_new
 
    p_prime = 0.5 * (Ap - Am)   # p' (정규화 단위)
    return p_prime
 
p_reflect = linearized_acoustic_pulse(bc='wall')
p_nscbc   = linearized_acoustic_pulse(bc='nscbc', sigma=0.05)
 
print(f"max |p'| (reflect): {np.max(np.abs(p_reflect)):.3e}")
print(f"max |p'| (nscbc):   {np.max(np.abs(p_nscbc)):.3e}")
# reflect: 0.998  /  nscbc: 0.025  → 두 자릿수 차이

같은 펄스에 대해 zero-gradient는 펄스를 거의 보존(반사)하고, NSCBC는 두 자릿수로 감쇠시킨다. $\sigma$ 만 바꿔 봐도 반사 계수가 거의 선형으로 따라온다.

실제 코드에 옮기기 전 체크리스트#

경계 종류부터 세라. subsonic outflow → BC 1개, supersonic outflow → BC 0개, subsonic inflow → BC 4개, supersonic inflow → BC 5개. 위 시각화로 직접 확인하자.
$\sigma$ 튜닝은 도메인 길이에 묶여 있다. 도메인을 2배 늘리면 $K$ 도 그에 따라 조정. 절대값으로 박지 말 것.
transverse 항을 그냥 떨어뜨리지 말 것. 3D LES에서는 Yoo & Im(2007)의 "relaxed" 항이 더 안정적이다.
충격이 출구에 닿는 case면 NSCBC만으로는 부족하다. Buffer zone(sponge layer)을 함께 써서 비선형 반사를 죽인다.
점성 플럭스는 별도로 처리한다. NSCBC가 자동으로 처리해 주지 않는다.

한 줄 정리#

출구는 "값"이 아니라 "들어오는 파의 진폭"을 정한다. 다섯 신호 중 들어오는 것 하나에 $\mathcal{L}_1 = K(p - p_\infty)$ 를 끼우면, 음파는 빠져나가고 평균은 표류하지 않는다. $\sigma$ 만 한 자릿수 안에서 튜닝하면, 6시간짜리 LES는 출구 때문에는 죽지 않는다.