초음속 유동은 상류를 모른다 — Euler 방정식의 특성곡선과 음파

영향권 쐐기가 수직선을 넘는 순간이 마하 1이고, 그 위에서는 출구에 걸 조건의 수가 특성 부호와 함께 줄어든다.

초음속 전투기는 자기 엔진 소리보다 먼저 도착한다. 머리 위를 지나간 뒤에야 굉음이 따라온다. 이유는 Euler 방정식의 특성속도에 있다. 작은 교란은 세 갈래 신호로 갈라지고, 그 자취가 $(x,t)$ 평면의 특성곡선(신호가 지나가는 자취)이 된다.

M 슬라이더를 0.4에 두면 파란 왼쪽 음파가 source의 상류로 뻗는다. 1을 넘기면 그 자취가 화면에서 사라진다.

background flow u₀ (Mach, c₀=1) = 0.40

The blue tint marks the upstream region the left-going sound has reached. Push u₀ past 1 and the blue tint vanishes: the C− signal now has a positive speed, so every signal runs downstream and the source can no longer talk to what lies upstream.

파란 선의 기울기 부호가 바뀌는 순간이다. 세 신호가 모두 오른쪽으로만 가는지, 왼쪽이 비는지 보면 된다.

압축성 방정식은 신호의 전파다#

1차원 비점성 압축성 유동은 세 보존식으로 적힌다.

\partial_t \rho + \partial_x(\rho u) = 0

\partial_t(\rho u) + \partial_x(\rho u^2 + P) = 0

\partial_t(\rho e_{\text{tot}}) + \partial_x\big[(\rho e_{\text{tot}} + P)\,u\big] = 0

$\rho$ 는 밀도, $u$ 는 속도, $P$ 는 압력, $e_{\text{tot}}$ 는 전체 비에너지다. 겉보기엔 유체가 "흐르는" 방정식처럼 보인다. 하지만 이 방정식의 본질은 이동이 아니라 신호의 전파다. 유동의 모든 동역학은 정해진 속도로 달리는 몇 개의 신호로 분해된다. 그 속도가 무엇인지 보려면 방정식을 잠깐 선형화하면 된다.

작은 교란이 파동 방정식을 만든다#

배경 상태 $(\rho_0, u_0, P_0)$ 위에 아주 작은 교란을 얹는다. $\rho = \rho_0 + \rho_1$ , $u = u_0 + u_1$ 로 두고 $P = K\rho^\gamma$ 을 쓴 뒤 1차 항만 남긴다. 배경을 따라 움직이는 미분 $D_t \equiv \partial_t + u_0\,\partial_x$ (공동이동 미분)를 정의하면, 두 교란식을 합쳐 하나의 파동 방정식이 나온다.

D_t^2\,\rho_1 - c_0^2\,\partial_x^2\,\rho_1 = 0

여기서 $c_0$ 는 음속이다.

c_0 = \sqrt{\gamma\,\frac{P_0}{\rho_0}}

$\gamma$ 는 비열비다. 이 방정식에 $\rho_1 \sim e^{i(kx-\omega t)}$ 를 넣으면 분산 관계가 떨어진다.

\frac{\omega}{k} = u_0 \pm c_0

$\omega$ 는 각진동수, $k$ 는 파수다. 교란은 배경 흐름 $u_0$ 에 실려서, 음속 $c_0$ 만큼 좌우로 퍼진다. 즉 한 점의 교란은 $u_0 + c_0$ 과 $u_0 - c_0$ 이라는 두 속도로 갈라진다.

세 갈래로 갈라지는 특성곡선#

두 음파가 전부가 아니다. 세 번째 신호는 유체 자신이다. 유체에 색을 칠했다고 하자. $t=0$ 에 왼쪽을 파랑, 오른쪽을 빨강으로 물들이면, 그 색 경계 $\phi$ 는 단순히 흐름을 따라 실려 간다.

\partial_t \phi + u_0\,\partial_x \phi = 0

이 신호는 속도 $u_0$ 으로 이동한다. 이것이 엔트로피(또는 접촉면) 특성이다. 1차원 유동의 특성 속도는 이렇게 세 개다.

\lambda_- = u_0 - c_0, \qquad \lambda_0 = u_0, \qquad \lambda_+ = u_0 + c_0

$(x,t)$ 평면에서 $\dfrac{dx}{dt} = \lambda$ 를 만족하는 선을 특성곡선이라 부른다. $\lambda_\pm$ 은 음파, $\lambda_0$ 은 물질의 이동이다. 각 음파 특성을 따라서는 Riemann 불변량 $u \pm \dfrac{2c}{\gamma-1}$ 이 일정하게 보존되고, 물질 특성을 따라서는 엔트로피가 보존된다.

$M = u_0/c_0$ 이 0.4이면 $\lambda_-$ 는 음수다. 왼쪽 음파는 상류로 간다. $M$ 이 1을 넘으면 $\lambda_- = u_0 - c_0$ 이 양수가 되고, 세 신호는 하류로만 달린다.

세 속도의 고유값 분해는 임의 방향 법선에서의 오일러 고유 구조로 가면 중복 고유값과 접선 성분이 붙는다. 1차원에서는 그 중복이 엔트로피 특성 하나로 접혀 있다.

초음속에서 신호의 원뿔이 넘어간다#

같은 이야기를 시공간 도표로 보면 더 선명하다. source에서 뻗어 나온 세 특성곡선 $C_-,\,C_0,\,C_+$ 사이의 쐐기는 영향권(domain of influence)이다. source가 신호를 보낼 수 있는 시공간 사건의 집합이다.

$M<1$ 이면 이 쐐기는 수직선을 사이에 두고 양쪽으로 벌어진다. 상류·하류 어디로든 신호가 간다. $M>1$ 이면 쐐기 전체가 오른쪽으로 넘어간다. $C_-$ 특성마저 오른쪽으로 기울어, 상류의 어떤 사건도 영향권 밖이다. 이것이 초음속 유동에서 "상류가 하류를 모르는" 이유다. 수치적으로도 초음속 출구에 경계조건을 걸면 안 되는 근거가 여기 있다.

아래 도표에서 둥근 탐침을 끌어 옮겨보자.

background flow u₀ (Mach, c₀=1) = 0.40

Drag the probe across the plane. At M = 0.4 the wedge straddles the vertical, so events both upstream and downstream are reachable. Slide u₀ above 1 and the whole wedge leans right — drop the probe anywhere to the left of the source and it stays red.

탐침을 쐐기 안에 놓으면 초록(도달 가능), 밖에 놓으면 빨강(도달 불가)으로 바뀐다. $u_0$ 을 1 위로 올린 뒤 탐침을 source 왼쪽 어디에 두어도 계속 빨강으로 남는다. 상류가 신호로부터 완전히 단절된 것이다.

Python으로 확인하는 특성 속도#

특성 속도가 정말 $u \pm c$ 와 $u$ 인지 확인해보자. 1차원 Euler를 원시 변수 $W=(\rho,u,P)$ 로 쓰면 준선형 형태 $\partial_t W + A(W)\,\partial_x W = 0$ 가 되고, 계수 행렬 $A$ 의 고유값이 곧 특성 속도다.

import numpy as np
 
def characteristic_speeds(rho, u, p, gamma=1.4):
    """1D Euler 원시변수 야코비안 A(W), W=(rho,u,p).
    A의 고유값 = 세 특성 속도."""
    c = np.sqrt(gamma * p / rho)              # 음속
    A = np.array([[u,      rho,        0.0  ],
                  [0.0,    u,          1/rho],
                  [0.0,    rho*c**2,   u    ]])
    lam = np.sort(np.linalg.eigvals(A).real)
    return c, lam
 
rho, u, p = 1.225, 220.0, 1.0e5               # 공기, u = 220 m/s
c, lam = characteristic_speeds(rho, u, p)
print(f"c0         = {c:8.2f} m/s")
print(f"eig(A)     = {np.round(lam, 2)}")
print(f"u-c,u,u+c  = {np.round([u-c, u, u+c], 2)}")
print(f"Mach       = {u/c:5.3f}  -> {'supersonic' if u > c else 'subsonic'}")

출력은 다음과 같다.

c0         =   338.06 m/s
eig(A)     = [-118.06  220.    558.06]
u-c,u,u+c  = [-118.06  220.    558.06]
Mach       =  0.651  -> subsonic

고유값이 정확히 $u-c,\,u,\,u+c$ 와 일치한다. 이번엔 실제로 교란을 시간전진시켜 신호의 앞면 속도를 재본다. 음향 Riemann 불변량 $w_\pm = u_1 \pm (c_0/\rho_0)\,\rho_1$ 은 각각 $u_0 \pm c_0$ 으로 실려 가는 스칼라라, 업윈드로 따로 옮기면 된다.

N, L = 400, 1.0
x = np.linspace(0, L, N, endpoint=False)
dx = L / N
rho0, c0 = 1.0, 1.0
u0 = 0.6 * c0                                  # 배경 Mach 0.6
 
bump = np.exp(-((x - 0.5) ** 2) / (2 * 0.01 ** 2))
wp = bump.copy()                              # 오른쪽 불변량 (u0 + c0)
wm = bump.copy()                              # 왼쪽  불변량 (u0 - c0)
 
def upwind_step(w, a, dt):
    if a > 0:
        dwdx = (w - np.roll(w, 1)) / dx
    else:
        dwdx = (np.roll(w, -1) - w) / dx
    return w - a * dt * dwdx
 
CFL = 0.4
dt = CFL * dx / (abs(u0) + c0)
steps = int(0.25 / dt)
for _ in range(steps):                        # 시간전진 루프
    wp = upwind_step(wp, u0 + c0, dt)
    wm = upwind_step(wm, u0 - c0, dt)
 
T = steps * dt
def front_speed(w):
    return (x[np.argmax(w)] - 0.5) / T
print(f"u0+c0: target {u0 + c0:+.3f}, measured {front_speed(wp):+.3f}")
print(f"u0-c0: target {u0 - c0:+.3f}, measured {front_speed(wm):+.3f}")

u0+c0: target +1.600, measured +1.600
u0-c0: target -0.400, measured -0.400

두 앞면이 각각 $u_0+c_0$ 과 $u_0-c_0$ 으로 움직인다. 오른쪽 불변량은 하류로, 왼쪽 불변량은 상류로. $u_0$ 을 $c_0$ 위로 올리면 두 번째 줄의 측정값도 양수로 바뀐다.

공기 숫자로 보면 간격이 더 선명하다. $\rho=1.225$ , $u=220$ m/s, $P=10^5$ , $\gamma=1.4$ 에서 음속은 338.06 m/s다. 출력의 고유값은 $-118.06$ , $220$ , $558.06$ 으로, $u-c$ , $u$ , $u+c$ 와 한 자리도 다르지 않다. $558.06-(-118.06)$ 은 676.12 m/s이고, 음속 338.06의 두 배와 같다. 가운데 220은 물질을 나르는 속도라 음속과 무관하다. 마하 0.651은 1보다 작다. 그래서 $\lambda_-$ 가 음수이고, 상류로 가는 신호가 하나 남아 있다.

두 번째 계산은 음속을 1로 둔 경우다. $u_0=0.6$ 이면 목표값은 $+1.600$ 과 $-0.400$ 이다. 측정값도 그 숫자다. $-0.400$ 은 $0.6-1$ 이다. 부호가 음수라는 점이 상류 도달을 뜻한다. $u_0$ 이 $c_0$ 보다 커지면 이 줄이 양수로 바뀌고, 시공간 쐐기가 오른쪽으로 넘어진 상태와 같아진다.

Riemann 불변량 $u \pm \dfrac{2c}{\gamma-1}$ 에서 $\gamma=1.4$ 이면 $\gamma-1=0.4$ 이고, $2/0.4=5$ 다. 음파 특성을 따라 일정한 양은 $u\pm 5c$ 다. 코드가 나르는 $w_\pm = u_1 \pm (c_0/\rho_0)\,\rho_1$ 은 이 양을 작은 교란으로 줄인 꼴이다. $\rho_0=c_0=1$ 이라 계수가 1로 남는다.

초음속 출구에 압력을 고정하면 파동이 되돌아온다#

영향권 쐐기가 수직선을 넘는 순간이 $M=1$ 이다. 그 위에서는 상류로 가는 특성이 없다. 출구가 초음속이면 세 특성의 부호가 모두 바깥쪽이다. 경계에서 정할 수 있는 값은 없다. 압력을 고정하면, 들어오지 않는 특성을 수치로 만들어 반사파가 도메인 안으로 들어온다.

아음속 출구는 다르다. $\lambda_-$ 만 안쪽을 향한다. 조건은 하나면 된다. 출구에서 되돌아오는 파동을 σ로 줄이는 NSCBC는 그 들어오는 특성에 반사 계수를 곱한다. 부호를 세기 전에 압력을 걸면, 마하 0.651의 공기에서도 상류 신호가 출구 처리와 섞인다.