출구는 열려 있는데 파동이 되돌아온다

비반사 출구에서 σ 하나가 반사율과 압력 표류를 동시에 정한다

출구는 열려 있다. 그런데 파동은 되돌아온다. 압축성 코드에서 출구 경계를 대충 외삽으로 처리하고 화염이나 난류를 계산하면, 도메인 한가운데에 이유 없는 진동이 자란다. 진동의 주기를 재보면 대개 영역 길이를 음속으로 나눈 값과 맞는다. 계산 영역 자체가 공명관이 되어 있다는 뜻이다. 오늘은 그 공명을 만드는 경계에서 무엇이 계산 가능하고 무엇이 지어낸 값인지, 그리고 그 지어낸 값을 조절하는 계수 하나가 어떤 대가를 요구하는지를 본다.

계산 영역의 끝은 물리가 아니다#

물리적 영역에는 끝이 없다. 연소기 출구 뒤에도 공간은 계속된다. 그런데 격자는 어딘가에서 끝나야 하고, 그 끝단 셀은 이웃이 없다. 이웃이 없으면 미분을 못 하고, 미분을 못 하면 지배방정식을 풀 수 없다.

RANS 코드에서는 이 문제가 오래 숨어 있었다. 난류 점성과 인공 점성이 크기 때문에, 경계에서 잘못 만들어진 파동은 몇 셀만 가도 소멸한다. LES와 DNS로 오면 사정이 달라진다. 인공 점성은 0에 가깝고 난류 점성도 최소다. 경계가 만든 오차는 사라지지 않고 영역을 가로질러 돌아다닌다.

Poinsot와 Lele이 1992년에 정리한 처방은 접근을 뒤집었다. 경계에서 변수를 외삽하는 대신, 경계를 가로지르는 파동을 세고 각각의 진폭을 결정한다. Euler 방정식에 대한 특성 경계조건(ECBC)을 점성항이 있는 Navier–Stokes로 확장한 것이어서 NSCBC(Navier–Stokes Characteristic Boundary Conditions)라고 부른다. 외삽 절차는 한 줄도 쓰지 않는다.

경계에서 셀 수 있는 것과 지어내야 하는 것#

경계가 $x_1 = L$ 에 있다고 하자. $x_1$ 방향 항들을 파동 형태로 다시 묶으면 연속방정식은 이렇게 된다.

\frac{\partial \rho}{\partial t} + d_1 + \frac{\partial (\rho u_2)}{\partial x_2} + \frac{\partial (\rho u_3)}{\partial x_3} = 0

$\rho$ 는 밀도, $u_i$ 는 속도 성분, $d_1$ 은 경계에 수직한 방향의 기여를 모은 항이다. 이 $d$ 벡터가 특성 해석의 산물이고, 그 안에는 파동 진폭 $\mathcal{L}_i$ 가 들어 있다.

\mathcal{L}_1 = (u_1 - c)\left( \frac{\partial p}{\partial x_1} - \rho c \frac{\partial u_1}{\partial x_1} \right)

\mathcal{L}_5 = (u_1 + c)\left( \frac{\partial p}{\partial x_1} + \rho c \frac{\partial u_1}{\partial x_1} \right)

$c$ 는 국소 음속( $c^2 = \gamma p / \rho$ ), $p$ 는 압력이다. $\mathcal{L}_1$ 은 $x_1$ 의 음의 방향으로 달리는 음파, $\mathcal{L}_5$ 는 양의 방향으로 달리는 음파의 진폭 변화율이다. 나머지 셋은 유체와 함께 실려 간다. $\mathcal{L}_2$ 는 엔트로피, $\mathcal{L}_3$ 과 $\mathcal{L}_4$ 는 접선 방향 속도 $u_2$ , $u_3$ 이고 모두 속도 $u_1$ 로 이동한다.

여기서 핵심은 속도의 부호다. 파동이 영역 밖으로 나가면 그 진폭은 내부 점들의 값으로 계산할 수 있다. 안으로 들어오면 그 정보는 해 안에 존재하지 않는다. 지어내야 한다. 들어오는 파동의 개수가 곧 그 경계에 줄 수 있는 물리적 경계조건의 개수다.

아래 다이어그램에서 Mach 수를 직접 밀어보자. 다섯 개의 특성선이 경계를 어느 방향으로 통과하는지가 실시간으로 바뀐다.

Mach M+0.30

boundary

$M$ 을 $+0.3$ 에서 $+1.4$ 로 올려보면 들어오던 $\mathcal{L}_1$ 이 방향을 바꾸고 필요한 조건 수가 1에서 0으로 떨어진다. 부호를 음수로 뒤집으면 같은 면이 유입구가 되면서 조건 수가 4로 뛴다. 아음속 출구에서 잘 돌던 코드가 초음속에서 발산한다면, 대개 이 표에서 허용하지 않는 조건을 하나 더 강제하고 있는 것이다.

LODI — 들어오는 파동의 진폭을 짓는 규칙#

들어오는 진폭을 지어내려면 근거가 필요하다. NSCBC는 경계의 각 점에서 접선 방향 항, 점성항, 반응항을 모두 지운 국소 1차원 비점성 시스템을 세운다. 이것이 LODI(Local One Dimensional Inviscid) 관계다.

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

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

LODI 관계는 물리적 조건이 아니다. 실제 계산에 쓰이는 방정식도 아니다. 오직 들어오는 $\mathcal{L}_i$ 를 추정하는 데만 쓴다. 절차는 세 단계다. 물리 조건이 걸린 보존방정식을 시스템에서 지우고, 지운 방정식에 대응하는 LODI 관계로 미지의 $\mathcal{L}_i$ 를 알려진 $\mathcal{L}_i$ 로 표현하고, 남은 방정식으로 나머지 변수를 시간 전진시킨다.

완전 비반사 출구는 여기서 가장 단순한 선택을 한다. 밖에서 들어오는 음파가 없다고 선언하는 것이다.

\mathcal{L}_1 = 0

들어오는 파동이 0이면 반사도 0이다. 깔끔해 보인다. 그런데 이 식 어디에도 바깥 압력 $p_\infty$ 가 없다.

σ = 0의 대가: 압력이 p∞로 돌아오지 않는다#

$p_\infty$ 가 없다는 것은 경계가 "지금 압력이 얼마여야 하는지"를 모른다는 뜻이다. 영역 안에서 열이 발생해 압력이 올라가도, 그 오프셋을 되돌릴 복원력이 어디에도 없다. 문제가 잘 정의되지 않는다.

Rudy와 Strikwerda의 처방은 들어오는 파동을 완전히 0으로 두는 대신, 압력 차이에 비례하도록 묶는 것이다.

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

$L$ 은 영역의 대표 길이, $M$ 은 최대 Mach 수, $\sigma$ 는 이 처방의 유일한 자유 파라미터다. $\sigma = 0$ 이면 완전 비반사로 되돌아간다. $\sigma$ 를 키우면 경계가 압력을 $p_\infty$ 쪽으로 끌어당긴다.

아래 시뮬레이션에서 직접 조작해보자. 왼쪽이 막힌 관, 오른쪽이 출구다. fire pulse로 압력파를 쏘고 sigma 슬라이더를 움직이면 된다.

sigma0.00

source q0.30

fire a pulse and watch the red curve at the outlet — that is the reflection. then drop sigma to 0 and watch the white curve settle above p_inf instead of returning to it.

볼 것은 두 가지다. 첫째, 펄스가 출구에 닿는 순간 붉은 곡선(들어오는 파동 $A^-$ )이 솟는지 — 이것이 반사다. 둘째, source q를 올린 채로 아래쪽 압력 이력을 보면, $\sigma = 0$ 에서는 흰 곡선이 $p_\infty$ 선으로 돌아오지 않고 위쪽에 머문다. 반사는 없앴지만 압력을 놓친 상태다.

두 실패 사이의 좁은 창#

$\sigma$ 는 서로 반대 방향으로 당기는 두 개의 실패 사이에 있다.

처방	출구에서 하는 일	실패 방식
B1 (외삽 + Riemann 불변량)	속도·밀도 외삽, 압력만 이완	외삽이 만든 가짜 파동
B2 (NSCBC, $\sigma = 0$ )	$\mathcal{L}_1 = 0$	평균 압력이 $p_\infty$ 에 묶이지 않음
B3 (NSCBC, $\sigma > 0$ )	$\mathcal{L}_1 = K(p - p_\infty)$	$\sigma$ 가 크면 반사
B4 (반사 출구)	압력 고정, $\mathcal{L}_1 = \mathcal{L}_5$	완전 반사 — 영역이 공명관이 됨

$\sigma$ 가 작으면 평균 압력이 떠내려가고, 크면 경계가 딱딱해져 음향 에너지를 되던진다. Poinsot와 Lele이 실제로 쓴 값은 $\sigma \approx 0.25$ , 외삽 기반 B1의 대응 계수로는 $\sigma = 0.58$ 이었다. 이 값들은 이론에서 유도된 것이 아니라 두 실패 사이에서 골라낸 타협점이다.

주파수 의존성을 보면 왜 타협인지 분명해진다. 각진동수 $\omega$ 인 음파에 대해 이 경계의 반사계수는

|R| = \frac{K}{\sqrt{K^2 + 4\omega^2}}

가 된다. 저주파일수록 반사가 크다. 즉 $\sigma$ 는 "낮은 주파수는 붙잡고 높은 주파수는 통과시키는" 필터다. 평균 압력은 $\omega \to 0$ 인 성분이므로 붙잡히고, 우리가 내보내고 싶은 음파는 통과한다. 이 분리가 통하는 구간이 좁은 창이다.

코드로 재보는 반사율과 압력 오프셋#

1차원 선형 음향이면 상태가 두 개의 특성 진폭으로 정확히 쪼개진다. $A^\pm = p' \pm \rho c u'$ 가 각각 $\pm c$ 로 이동한다. $\Delta t = \Delta x / c$ 로 잡으면 매 스텝 정확히 한 칸씩 밀면 되므로, 화면에 보이는 흔들림은 전부 경계조건이 만든 것이다.

import numpy as np
 
N, C, L = 240, 1.0, 1.0
DX = L / N
DT = DX / C                       # 정확히 한 칸 이동 — 스킴 확산 없음
 
 
def outlet_relax_k(sigma, mach=0.0):
    """NSCBC 이완계수 K = sigma (1 - M^2) c / L"""
    return sigma * (1.0 - mach ** 2) * C / L
 
 
def duct_step(ap, am, am_b, k_relax, q):
    """A+ 는 오른쪽 한 칸, A- 는 왼쪽 한 칸. 출구에서 A- 는 지어내야 한다."""
    ap[1:] = ap[:-1].copy()
    ap[0] = 0.0
    am[:-1] = am[1:].copy()
 
    p_b = 0.5 * (ap[-1] + am[-1])
    am_b -= DT * k_relax * p_b     # L1 = K (p - p_inf)
    am[-1] = am_b
 
    ap[0] = am[0]                  # 왼쪽은 막힌 끝 (u = 0)
    ap += q * DT                   # 약한 균일 열발생
    am += q * DT
    return am_b
 
 
def measure_outlet(sigma, q, steps, pulse):
    ap, am, am_b = np.zeros(N), np.zeros(N), 0.0
    x = (np.arange(N) + 0.5) * DX
    if pulse:
        ap += np.exp(-((x - 0.30) / 0.09) ** 2)
    k = outlet_relax_k(sigma)
    refl = 0.0
    for n in range(steps):
        am_b = duct_step(ap, am, am_b, k, q)
        if pulse and n > 0.85 * N:
            refl = max(refl, np.abs(am[:-8]).max())
    return refl, float(np.mean(0.5 * (ap + am)))
 
 
for sigma in (0.0, 0.25, 1.0, 4.0, 10.0):
    r, _ = measure_outlet(sigma, q=0.0, steps=650, pulse=True)
    _, p = measure_outlet(sigma, q=0.3, steps=6000, pulse=False)
    print(f"sigma={sigma:5.2f}   반사={r * 100:5.1f}%   평균압력-p_inf={p:+.4f}")

실행 결과다.

sigma= 0.00   반사=  0.0%   평균압력-p_inf=+0.3000
sigma= 0.25   반사=  1.9%   평균압력-p_inf=+0.0007
sigma= 1.00   반사=  7.3%   평균압력-p_inf=+0.0006
sigma= 4.00   반사= 24.5%   평균압력-p_inf=+0.0227
sigma=10.00   반사= 46.8%   평균압력-p_inf=-0.1142

$\sigma = 0$ 은 반사를 완전히 없애지만 열발생이 만든 오프셋 0.3을 그대로 남긴다. $\sigma = 0.25$ 에서 반사는 2% 미만이고 압력은 $p_\infty$ 에서 0.0007 이내로 잡힌다. 여기까지는 예상대로다.

예상 밖은 마지막 두 줄이다. $\sigma$ 를 4, 10으로 키우면 반사가 25%, 47%로 뛰는 것까지는 당연한데, 압력 오프셋도 다시 나빠진다. 경계가 너무 뻣뻣해지면 스스로 진동을 만들고, 그 진동이 평균값을 흔든다. $\sigma$ 를 키우는 것은 압력을 붙잡는 대가로 반사를 사는 거래가 아니다. 어느 지점을 넘으면 둘 다 잃는다.

격자가 만든 파동은 반대 방향으로 되돌아온다#

경계에서 반사되는 것이 물리적 음파만은 아니다. 파장이 격자 간격의 네 배보다 짧은 성분은 물리적 해가 아니라 이산화가 만든 것이다. Poinsot와 Lele은 이것을 "q 파동"이라 불러 물리적인 "p 파동"과 구분했다.

q 파동의 특징은 군속도다. 이동 속도 $V$ 인 1차원 이류 방정식에서도 짧은 파장의 군속도 $u_g$ 는 $V$ 와 부호가 반대다. 유동이 오른쪽으로 흐르는데 수치 오차는 왼쪽으로 거슬러 올라간다. 게다가 $|u_g / V|$ 는 스킴의 차수가 높을수록 커진다. 고차 스킴이 더 위험하다는 뜻이다.

그래서 경계조건은 두 개의 반사계수로 평가해야 한다. 물리 파동의 반사 $A_p / A_1$ 과 수치 파동의 반사 $A_q / A_1$ 이다. 쓸 만한 경계조건이라면 어떤 경우에도 $A_q / A_1 \ll 1$ 이어야 하고, 비반사를 표방한다면 $A_p / A_1 \ll 1$ 까지 만족해야 한다. 계산 시작 시 초기장에 가파른 구배를 넣는 것만으로도 q 파동은 생성되고, DNS에서는 소멸하지 않은 채 남는다.

출구 조건을 고를 때의 한 줄#

$\sigma$ 는 튜닝 노브가 아니라 두 개의 실패 사이에 놓인 좌표다. 왼쪽 끝에는 어디에도 묶이지 않은 압력이, 오른쪽 끝에는 공명관이 된 계산 영역이 있다. 0.25 근처가 권장되는 이유는 그 지점에서 저주파만 붙잡고 나머지를 통과시키기 때문이다.

출구에서 원인 모를 진동을 만났다면 순서는 이렇다. 먼저 그 경계에서 들어오는 특성 파동의 개수를 세고, 지금 강제하는 조건 수가 그 값과 같은지 확인한다. 다음으로 진동 주기가 $2L/c$ 의 배수인지 본다. 맞으면 경계 반사이지 물리가 아니다. 마지막으로 $\sigma$ 를 낮춰본다. 진동이 줄면 원인은 경계였고, 평균 압력이 떠내려가기 시작하면 반대편 실패로 넘어간 것이다.

출구는 열려 있는데 파동이 되돌아온다 — NSCBC와 σ가 정하는 반사