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

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

Mach 슬라이더를 양수에서 음속으로 넘기면 들어오던 화살표가 나간다. 부호를 음수로 두면 같은 면의 조건 개수가 뛴다.

아음속 출구의 조건은 하나다. 초음속으로 넘기면 0, 유입으로 뒤집으면 4가 된다. 화살표가 영역 안을 향하는 동안만 그 파동은 지어낸 값이다.

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

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

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

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

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

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

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

$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$로 이동한다.

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

특성선이 음파를 영역 전체에서 어떻게 나누는지는 초음속 유동이 상류를 모르는 자리에 있다. 여기서는 그 선이 출구 한 면에서 끊기는 개수만 센다.

$M$을 $+0.3$에서 $+1.4$로 올리면 앞에서 뒤집힌 화살표가 $\mathcal{L}_1$이다. 조건 수는 1에서 0. 부호를 음수로 두면 유입구가 되고 조건은 4다. 아음속 출구에서 잘 돌던 코드가 초음속에서 발산한다면, 이 개수를 넘는 조건을 하나 더 걸고 있는 경우가 많다.

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

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

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

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

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

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

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

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

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

막힌 관의 sigma와 source q는 출력 표 다음에 둔다.

두 실패 사이의 좁은 창#

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

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

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

가 된다. 저주파일수록 반사가 크다. 즉 $\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

fire pulse로 쏜 펄스가 출구에 닿을 때 붉은 곡선 $A^-$가 솟으면 반사다. source q를 올린 채 $\sigma = 0$이면 흰 곡선이 $p_\infty$ 위에 머문다. 반사는 없고, 압력은 놓친 상태다.

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

예상 밖은 마지막 두 줄이다. $\sigma$를 4, 10으로 키우면 반사가 25%, 47%로 뛰는 것까지는 당연한데, 압력 오프셋도 다시 나빠진다. 경계가 너무 뻣뻣해지면 스스로 진동을 만들고, 그 진동이 평균값을 흔든다. $\sigma$를 키우는 것은 압력을 붙잡는 대가로 반사를 사는 거래가 아니다. 어느 지점을 넘으면 둘 다 잃는다. 출력의 반사 열은 0.0, 1.9, 7.3, 24.5, 46.8%다. 코드는 $C = 1$, $L = 1$, Mach 0이라 $K = \sigma \left( 1 - M^2 \right) \frac{c}{L}$이 $K = \sigma$로 줄어든다. σ를 0.25에서 10으로, 40배 키울 때 반사는 1.9에서 46.8이다. 약 25배. 선형이 아니다. $|R| = \frac{K}{\sqrt{K^2 + 4\omega^2}}$는 K가 커질수록 1로 포화한다. 펄스는 단일 주파수가 아니라 폭 0.09의 가우스다. 한 ω로 이 퍼센트를 역산하지는 않는다. 방향만 맞으면 된다. K를 키우면 반사가 늘고, 늘다가 둔해진다.

압력 열은 그 단조를 따르지 않는다. $q = 0.3$인 실행에서 $\sigma = 0$의 평균은 +0.3000, 열원과 같은 크기다. $\sigma = 0.25$에서 +0.0007, $\sigma = 1$에서 +0.0006. 0.3000 대비 0.0007은 0.23%다. 여기까지 이완이 오프셋을 거의 지운다. $\sigma = 4$에서 +0.0227. 0.0007의 약 32배다. $\sigma = 10$에서 −0.1142. 부호가 뒤집힌다. 0.0007 대비 절대값은 약 160배다. 반사 46.8%인 경계는 평균을 붙잡지 못하고 흔든다. 표에서 오프셋이 바닥인 구간은 0.25와 1.0이다. Poinsot와 Lele이 쓴 0.25는 그 바닥의 왼쪽 끝이다. 0.58은 B1 대응값이라 이 열에 그대로 찍히지 않는다.

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

경계에서 반사되는 것이 물리적 음파만은 아니다. 파장이 격자 간격의 네 배보다 짧은 성분은 물리적 해가 아니라 이산화가 만든 것이다. 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에서는 소멸하지 않은 채 남는다.

주기가 2L/c면 출구의 σ를 먼저 의심한다#

$\sigma$는 튜닝 손잡이가 아니다. 한쪽 끝은 어디에도 안 묶인 압력, 다른 끝은 공명관이 된 영역이다. 0.25 근처는 저주파만 붙잡고 나머지를 통과시키는 자리다.

원인 모를 진동이면 들어오는 특성 파동의 개수와 강제 조건 수가 같은지 센다. 주기가 $2L/c$의 배수면 물리가 아니라 경계다. $\sigma$를 낮춰 진동이 줄면 경계가 원인이고, 평균 압력이 떠내려가면 반대편으로 넘어간 것이다. 같은 출구를 파동 개수부터 세는 쪽은 NSCBC를 비반사 출구로 세운 이전 정리에 있다.