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

충격 속도는 두 상태의 산술평균이다. 왼쪽이 더 작으면 그 점프는 충격이 아니라 팽창파다.

1948년 Burgers는 난류 모형으로 방정식 한 줄을 던졌다. 매끈한 곡선이 외부 힘 없이 유한 시간 안에 절벽이 된다. 특성선(characteristic, 정보가 흐르는 시공간 곡선)이 겹치는 시각이 그 절벽의 생일이다. 약해(weak solution, 미분이 깨진 뒤 적분형으로 남는 해)는 그 겹침을 한 줄로 자른다.

진폭 $A$ 와 파수 $k$ 를 올리면 빨간 점선, 깨짐 시각이 위로 당겨진다.

Amplitude A: 0.35Wave number k: 1

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

봉우리 직선이 앞 골짜기 직선을 덮는 순간이 그 점선이다. $A$ 가 커지면 만남이 앞당겨지고, $k = 2$ 면 두 곳에서 동시에 겹친다.

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

쌍곡 방정식의 본질은 신호 전파다. 정보가 특성선을 따라 일정한 속도로 흐르는 한 모든 게 평화롭다. 선형 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^*$ 를 넘는 순간 수직 절벽이 솟는다.

기울기 최솟값은 $-2\pi k A$ 다. $2\pi \approx 6.28$ 이니 $2\pi \times 0.35 \approx 2.20$ 이고, 역수가 바로 $0.455$ 다. 봉우리는 $0.55 + 0.35 = 0.90$ , 골짜기는 $0.20$ 이다. 처음 높이 $0.70$ 이 나중 shock height 약 $0.7$ 과 맞는다. $t = 1.2$ 까지 가도 점프가 새로 자라지는 않는다. 모서리만 둥글어진다.

$A$ 를 $0.45$ 로 올리면 분모가 약 $2.83$ 이 되고 $t^*$ 는 $0.35$ 근처로 준다. 파수까지 $2$ 로 두면 분모가 두 배, $t^* \approx 0.227$ 이다. 진폭이 같아도 파장이 짧으면 더 빨리 겹친다.

약해와 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이 답으로 살아남는다.

$s = (q_L + q_R)/2$ 에 $0.90$ 과 $0.20$ 을 넣으면 $s = 0.55$ 다. 사인파의 배경값과 같다. 충격은 그 평균 속도로 미끄러지고, 앞뒤 상태만 남는다. 순서를 뒤집어도 속도 공식은 $0.55$ 를 내놓는다. 답만 충격에서 팽창파로 바뀐다. 같은 산술평균이 두 파를 구분하지 못한다.

격자 간격은 $1/256 \approx 0.00391$ 이다. CFL $0.4$ 와 초기의 $\max|q| = 0.90$ 이면 $\Delta t \approx 0.00174$ 다. $t = 1.2$ 까지 약 $690$ 스텝이다. 충격 뒤에서 $\max|q|$ 가 국소적으로 커지면 같은 식이 걸음을 더 잘게 쪼갠다. 셀 경계의 그 칼은 충격관 한 칸에서 다시 만난다.

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 문제의 정확해를 셀 인터페이스마다 계산해 약해를 보장한다.

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

앞에서 민 그림에서 청록 곡선과 흐린 초기 점선이 충격 좌우로 잘라낸 면적은 같다. equal area rule이다. $k = 2$ 에서 겹침이 두 곳이면, 두 충격은 곧 하나로 합쳐지며 점프가 커진다.

그 면적이 같다는 말은, 자른 뒤에도 $\int q\,dx$ 가 양끝 플럭스 차이와 맞는다는 뜻이다. 봉우리 $0.90$ 과 골짜기 $0.20$ 의 평균 $0.55$ 로 충격이 걸어가면, 잘린 삼각형 둘이 서로 상쇄된다. 높이 $0.70$ 은 둔 채 위치만 $s = 0.55$ 로 옮긴다.

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

(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 속도로 뒤로 밀린다.

정체를 연속 곡선으로 붙들고 있으면 다중값이 남는다. 한 줄로 자를지, 팽창으로 풀지는 왼쪽과 오른쪽의 순서가 가른다. 계단에서 진동이 보이면 플럭스 리미터가 먼저 의심할 자리다.