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

초음속 전투기는 자기 엔진 소리보다 먼저 도착한다. 머리 위를 지나간 뒤에야 굉음이 따라온다. 왜 이런 일이 벌어지는지는 사실 Euler 방정식 안에 이미 적혀 있다. 이 글은 작은 교란이 어떻게 세 갈래 신호로 갈라지는지 유도한다. 그 신호가 $(x,t)$ 평면에 그리는 특성곡선(신호가 지나가는 자취)을 조작 가능한 그림으로 확인하고, 마지막엔 세 특성 속도를 Python으로 직접 뽑아본다.

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

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}$ 이 일정하게 보존되고, 물질 특성을 따라서는 엔트로피가 보존된다.

아래 애니메이션에서 배경 흐름 $u_0$ 을 직접 조작해보자.

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.

$M = u_0/c_0$ 을 0.4에 두면 파란 신호(왼쪽 음파)가 source의 왼쪽, 즉 상류로 뻗어간다. $M$ 을 1 위로 올리면 그 파란 자취가 사라진다. $\lambda_- = u_0 - c_0$ 이 양수로 바뀌면서 세 신호가 모두 하류로만 달리기 때문이다.

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

같은 이야기를 시공간 도표로 보면 더 선명하다. 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$ 위로 올리면 두 번째 줄의 측정값도 양수로 바뀐다.

특성곡선이 남기는 것#

압축성 유동의 동역학은 세 특성 속도 $u-c,\ u,\ u+c$ 로 달리는 신호의 전파다. 두 개는 음파, 하나는 물질의 이동.
영향권 쐐기가 수직선을 넘어가는 순간이 $M=1$ 이다. 그 위에서는 상류로 가는 신호가 없다.
초음속 경계에서 특성이 모두 한 방향이면, 걸 수 있는 경계조건의 수도 그만큼 정해진다. 특성의 부호를 세는 습관이 경계조건 실수를 막는다.