특성선이 부딪히면 충격파가 태어난다 — 버거스 방정식과 (x,t) 평면

1948년, 네덜란드 물리학자 J. M. Burgers는 난류의 모형으로 단 한 줄짜리 방정식을 던졌다. 매끈한 곡선이 외부 힘 없이도, 유한한 시간 만에 수직으로 깎이며 충격파가 되는 — 비선형 쌍곡 PDE의 가장 짧은 드라마. 이 글은 특성선(characteristic, 정보가 흐르는 시공간 곡선)의 기하학으로 충격파의 탄생 시각을 손으로 계산하고, 약해(weak solution, 미분 가능성이 깨진 곳에서 적분형으로 정의되는 해)가 어떻게 다중값을 한 줄로 자르는지 정리한 뒤, Python 60줄로 같은 그림을 재현한다. 마지막엔 직접 직선이 부딪히는 (x,t) 평면을 만질 수 있다.

직선들이 같은 점에서 만난다#

쌍곡 방정식의 본질은 신호 전파다. 정보가 특성선을 따라 일정한 속도로 흐르는 한 모든 게 평화롭다. 선형 advection $\partial_t q + u\,\partial_x q = 0$ 에서 $u$ 가 상수면 특성선은 평행한 직선이고, 초기 모양이 그대로 옆으로 이동한다. 이 그림에는 충격파가 들어설 자리가 없다.

비선형성이 등장하는 순간 모든 게 달라진다. $u$ 가 $q$ 자체에 의존하면, 큰 값은 빠르게, 작은 값은 느리게 간다. 빠른 값이 느린 값을 추월하는 순간 (x,t) 평면 한 점에 두 개의 정보가 동시에 도착한다. 함수 값이 다중값(multi-valued)이 되는 순간, 즉 물리적으로 말이 안 되는 그 순간이 바로 충격파의 탄생이다.

자기 자신을 운반하는 방정식#

버거스 방정식의 보존형은 다음과 같다.

\partial_t q + \partial_x \left(\tfrac{1}{2} q^2\right) = 0

$q$ 는 보존되는 스칼라, $f(q) = q^2/2$ 는 플럭스다. 체인룰을 풀면 비보존형이 된다.

\partial_t q + q\,\partial_x q = 0

advection 속도 $u(x,t) = q(x,t)$ . 자기 자신이 자기 자신을 운반하는 구조다. 이 비선형성이 모든 드라마의 출발점이다. 자코비안 $\partial f/\partial q = q$ 가 양수인 곳에선 신호가 오른쪽으로, 음수인 곳에선 왼쪽으로 간다. 부호가 자리마다 다른 초기 분포는 — 사인 곡선처럼 — 빠른 곳이 느린 곳을 따라잡게 만든다.

특성선이 그리는 (x,t) 평면#

특성선은 $dx/dt = q$ 를 만족하는 곡선이다. 그 위에서 동행 미분 $Dq/Dt$ 가 0이므로, $q$ 는 특성선을 따라 출발 값으로 동결된다. 결국 특성선은 (x,t) 평면에서 기울기 $1/q_0(x_0)$ 의 직선이다. 출발점 $x_0$ 의 초기 값 $q_0(x_0)$ 가 그 직선의 속도가 된다.

초기 조건이 $q_0(x_0) = 0.55 + A\sin(2\pi k x_0)$ 일 때:

봉우리( $q_0 > 0.55$ )에서 출발한 직선은 빠르게 오른쪽으로 기운다
골짜기( $q_0 < 0.55$ )에서 출발한 직선은 느리게 기운다

봉우리 직선이 바로 앞 골짜기 직선의 발등을 밟는 시각이 반드시 온다. 두 직선이 한 점에서 만나는 그 점이 첫 충격이다.

깨지는 시각 t* — 처음 부딪히는 순간#

이웃한 출발점 $x_0$ , $x_0 + \Delta x$ 의 특성선이 시각 $t$ 에 만나려면 다음이 성립해야 한다.

x_0 + q_0(x_0)\,t = (x_0 + \Delta x) + q_0(x_0 + \Delta x)\,t

$\Delta x \to 0$ 극한을 취하면:

1 + q_0'(x_0)\,t = 0 \;\;\Longrightarrow\;\; t = -\frac{1}{q_0'(x_0)}

음의 기울기일 때만 만난다. 가장 이른 시각은 $q_0'(x_0)$ 가 가장 음수인 점에서 일어난다.

t^* = \frac{1}{\max_{x_0}\big(-q_0'(x_0)\big)}

$q_0$ 의 진폭과 주파수가 클수록 더 빨리 깨진다. 위 사인 함수의 경우 $q_0'(x_0) = 2\pi k A\cos(2\pi k x_0)$ , 최소값은 $-2\pi k A$ 이다.

t^* = \frac{1}{2\pi k A}

$A = 0.35$ , $k = 1$ 이면 $t^* \approx 0.455$ . 이 시각 이전엔 곡선이 부드럽게 비대칭화되고, $t^*$ 를 넘는 순간 수직 절벽이 솟는다.

약해와 Rankine–Hugoniot — 다중값을 자르는 칼#

$t > t^*$ 에서 직진 해석은 한 점에 세 개의 $q$ 를 동시에 토해낸다. 물리적으로 불가능하다. 해법은 약해다. 미분형 대신 적분 보존형으로 정의를 약화시킨다.

\frac{d}{dt}\int_{x_L}^{x_R} q\, dx \;=\; f\big(q(x_L)\big) - f\big(q(x_R)\big)

이 식은 $q$ 가 미분 불가능해도 — 즉 불연속이어도 — 의미가 있다. 두 평탄 상태 $q_L$ , $q_R$ 이 충격으로 연결되어 속도 $s$ 로 전파된다고 두면 다음이 따라온다.

s = \frac{f(q_L) - f(q_R)}{q_L - q_R} = \frac{q_L + q_R}{2}

이것이 Rankine–Hugoniot 조건이다. 충격 속도는 두 상태의 산술 평균이다. 이 조건이 다중값 영역을 한 줄로 자르고, 자른 양쪽 면적이 같도록 — 흔히 "equal area rule"이라 부른다 — 만든다.

추가로 entropy 조건이 필요하다. $q_L > q_R$ 일 때만 충격이 물리적이고, $q_L < q_R$ 이면 부채꼴 rarefaction(확장파)으로 풀어야 한다. 이 조건을 어기면 수학적으론 약해이지만 비물리적인 expansion shock이 답으로 살아남는다.

Python: 특성선 추적부터 Godunov 충격까지#

특성선을 직접 그리고, Godunov 1차 스킴으로 약해 충격을 잡아본다.

import numpy as np
 
def burgers_breaking_time(q0_prime_min):
    # q0_prime_min: x에 대한 q0'의 최솟값 (음수여야 충격 발생)
    if q0_prime_min >= 0:
        return float('inf')
    return -1.0 / q0_prime_min
 
def trace_characteristics(q0_func, x0_arr, t_arr):
    # x[i, j] = x0_arr[i] + q0(x0_arr[i]) * t_arr[j]
    return x0_arr[:, None] + q0_func(x0_arr)[:, None] * t_arr[None, :]
 
def rankine_hugoniot_speed(qL, qR):
    # 충격 속도 s = (qL + qR) / 2
    return 0.5 * (qL + qR)
 
def godunov_burgers_step(q, dx, dt):
    # 보존형: q^{n+1} = q^n - (dt/dx) * (F_{i+1/2} - F_{i-1/2})
    qL = q                          # 셀 좌측 = 현재 셀 값
    qR = np.roll(q, -1)             # 셀 우측 = 다음 셀 값
    s = 0.5 * (qL + qR)
    f_shock = np.where(s >= 0, 0.5 * qL**2, 0.5 * qR**2)
    f_rare  = np.where(qL >= 0, 0.5 * qL**2,
              np.where(qR <= 0, 0.5 * qR**2, 0.0))
    F_iph = np.where(qL > qR, f_shock, f_rare)   # 우측 인터페이스 플럭스
    F_imh = np.roll(F_iph, 1)                    # 좌측 인터페이스 = 한 칸 시프트
    return q - (dt / dx) * (F_iph - F_imh)
 
# 시뮬레이션 파라미터 -----------------------------------------
N, L = 256, 1.0
A, k = 0.35, 1
x = (np.arange(N) + 0.5) * (L / N)
q0 = lambda xx: 0.55 + A * np.sin(2 * np.pi * k * xx)
q  = q0(x)
 
t_star = burgers_breaking_time(-2 * np.pi * k * A)
print(f"breaking time t* = {t_star:.3f}")
 
t_end, t = 1.2, 0.0
while t < t_end:
    dt = 0.4 * (L / N) / np.max(np.abs(q))    # CFL = 0.4
    q  = godunov_burgers_step(q, L / N, dt)
    t += dt
print(f"shock height (max - min) at t = {t:.2f}: {q.max() - q.min():.3f}")

breaking time t* = 0.455, shock height ≈ 0.7이 출력된다. Godunov 스킴은 다중값 영역을 자동으로 한 줄로 잘라준다 — 그 칼이 바로 위 np.where(qL > qR, f_shock, f_rare) 한 줄이다. 이 한 줄이 Riemann 문제의 정확해를 셀 인터페이스마다 계산해 약해를 보장한다.

직접 직선이 부딪히는 걸 보자#

아래 시뮬레이션에서 진폭과 파수를 바꿔보자.

Amplitude A: 0.35Wave number k: 1

q0(x) = 0.55 + A sin(2 pi k x)

진폭 $A$ 를 0.45까지 올리면 $t^*$ 가 0.35 근처로 짧아지고, 빨간 점선이 위로 올라간다. 파수 $k$ 를 2로 하면 깨짐이 두 곳에서 동시에 일어나며, 일단 깨진 뒤엔 두 충격이 합쳐지면서 더 큰 충격으로 자란다. 위 패널의 청록 곡선과 흐린 점선(초기 곡선) 사이의 면적이 같다는 사실 — equal area rule — 을 직접 눈으로 확인할 수 있다.

수치 충격 캡처가 실패하는 두 가지 이유#

(1) CFL 위반. $\Delta t \le \Delta x / \max|q|$ 를 어기면 수치 영역이 물리 영역을 추월한다. 진폭이 시간에 따라 자라는 것 같으면 dt부터 의심하라. 위 코드처럼 매 스텝 np.max(np.abs(q))로 dt를 다시 계산하는 게 안전하다. 충격이 들어선 직후 $q$ 가 국소적으로 커질 수 있으므로 더더욱 그렇다.

(2) Entropy 잘못 잡기. Rankine–Hugoniot만 만족시키면 비물리적인 expansion shock도 답이 된다. Roe 같은 단순 선형화 스킴은 transonic rarefaction(부채꼴이 0을 가로지를 때) 영역에서 "entropy glitch"를 만든다. Godunov, HLLE, exact Riemann solver는 자동으로 entropy 조건을 보장하지만, Roe를 쓰려면 Harten 식 entropy fix가 명시적으로 필요하다.

LES 압축성 코드에서 충격 위치가 떨리거나 음의 압력이 보이면 거의 항상 위 둘 중 하나다.

다음에 정체가 시작되거든#

자동차들이 한 차로에 모여 속도가 떨어진다. 앞 차는 느리고, 뒤 차는 빠르다. 빠른 차가 느린 차의 범퍼에 닿는 순간 — 정체의 머리가 만들어진다. (x,t) 평면에서 자동차 궤적이 만나는 점이다. 도로 정체는 사실상 버거스 방정식의 충격파고, 풀리지 않는 한 Rankine–Hugoniot 속도로 뒤를 향해 전파한다.

한 줄 정리. 비선형 쌍곡 PDE는 부드러운 초기값에서도 유한 시간에 불연속을 만든다. 그 사실을 받아들이고, 약해라는 우회로를 통해 의미를 회복하는 것이 충격파 수치해석의 출발점이다.