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

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

Mach 슬라이더를 0에서 1.6까지 밀어 보자. 숫자가 오르면 안으로 들어오는 화살표가 줄어든다.

아음속에서는 굵은 빨간 화살표가 하나 남는다. 1을 넘기면 그 화살표가 경계를 밖으로 나간다. 출구가 지정할 신호의 개수가 그 화면에 있다.

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

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

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

출구가 보낸 다섯 신호#

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

여기서 $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 슬라이더에서 굵은 선으로 보이던 것이 $u-c$다. 아음속 출구는 그 하나만 안으로 들어온다. 필요한 경계조건은 하나다. 압력·속도·온도를 한꺼번에 걸면 over-specified가 되어 솔버가 죽는다.

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

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

여기서 $\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으로 둔다.

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

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

$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 사이에서 튜닝한다.

통과 시간 $\tau$를 도메인 횡단 $L/c$로 두면 $L/(c\tau)=1$이라 $|R|\approx\sigma$가 된다. 추천값 $\sigma\approx 0.25$는 이 추정에서 진폭 반사 약 0.25다. $\sigma$를 0.1에서 0.5로 옮기면 $|R|$도 그 구간을 따라간다. $K$는 $\sigma$에 비례한다. 0.25에서 0.5로 두 배 하면 $K$도 두 배다. $(1-M_{\max}^{2})$는 마하 보정이다. $M_{\max}=0$이면 인자는 1, $M_{\max}=0.5$이면 $1-0.25=0.75$라 $K$는 25% 작아진다. $M_{\max}$가 1에 붙으면 $K$는 0으로 꺼지고, 들어오는 파가 사라지는 그림과 같은 방향이다.

모드를 Zero-gradient와 NSCBC로 바꿔 보자. $\sigma$ 슬라이더는 들어오는 파에 곱하는 계수다.

Zero-gradient에서는 가우시안 펄스가 오른쪽 경계에서 거의 그대로 튕긴다. 반사 계수는 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$만 바꿔 봐도 반사 계수가 거의 선형으로 따라온다.

주석의 피크는 0.998과 0.025다. 비는 약 40이다. 코드의 $\sigma=0.05$는 나가는 파에 곱해지는데, 피크 비는 그 0.05보다 더 작은 0.025 쪽에 가깝다. 0.3 bar를 6시간으로 나누면 0.05 bar/h다. $\mathcal{L}_1=0$의 $\partial p/\partial t=-\mathcal{L}_5/2$가 그 기울기의 부호를 정한다.

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

• 경계 종류부터 세라. 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가 자동으로 처리해 주지 않는다.

σ를 0에 붙이면 평균이 걷는다#

$\mathcal{L}_1=0$은 음파를 깨끗이 내보내지만, $\partial p/\partial t=-\mathcal{L}_5/2$가 남는다. $\mathcal{L}_5$의 시간평균이 작은 양수면 압력은 상한이 없다. 6시간 뒤의 0.3 bar는 그 적분이 로그에 남긴 자국이다. $\sigma\approx 0.25$를 $K=\sigma(1-M_{\max}^{2})c/L$에 넣으면 표류를 당기되, $\tau=L/c$인 추정에서는 진폭의 약 1/4을 되돌린다.

$\sigma$를 5 근처로 올리면 경계가 다시 딱딱해져 펄스가 되튄다. 0과 5 사이, 코드가 쓴 0.05와 추천 0.25는 같은 직선의 두 점이다. 여기서는 들어오는 파가 하나라는 데까지 왔다. 그 한 개의 계수를 더 조이는 쪽은 출구의 σ가 반사를 어떻게 고정하는지에 있다.