충격파가 제자리에 멈춰 섰다 — 보존형과 원시형이 갈라지는 자리

미분 규칙으로는 같은 두 식이, 불연속에서는 다른 물리를 푼다. 검사체적에서 유도된 원래 형태를 유지한 쪽만 옳은 속도를 낸다.

1차원 Burgers 솔버를 두 벌 짜서 붙여 본 적이 있다. 한쪽은 flux 차분, 다른 쪽은 속도에 기울기를 곱한 형태다. 종이 위에서 두 식은 연쇄법칙 한 줄로 서로 옮겨 간다. 그런데 리만 문제를 넣자 한쪽 충격파가 아예 움직이지 않았다. 이 글은 그 정지의 원인을 검사체적 유도까지 되짚고, 격자를 8배 세분해도 회복되지 않는다는 것을 Python으로 확인한다.

충격파가 제자리에 멈춰 섰다#

같은 방정식을 두 가지로 쓸 수 있다. 보존형(conservative form)은 시간 변화율과 flux 발산으로 쓴다.

\frac{\partial u}{\partial t} + \frac{\partial}{\partial x}\left(\frac{u^2}{2}\right) = 0

원시형(primitive form, 비보존형)은 미분을 풀어 속도에 기울기를 곱한다.

\frac{\partial u}{\partial t} + u\,\frac{\partial u}{\partial x} = 0

$u$ 가 매끄러우면 $\partial_x(u^2/2) = u\,\partial_x u$ 이므로 두 식은 같다. 곱미분 규칙이 그렇게 말한다.

리만 문제를 넣어 보자. 왼쪽 $u_L = 1$ , 오른쪽 $u_R = 0$ 이다. 정확해는 속도 $0.5$ 로 오른쪽으로 가는 충격파다. 보존형 Godunov 스킴은 $0.5006$ 을 냈다. 원시형 풍상차분(upwind difference)은 $0.0000$ 을 냈다. 충격파가 출발선에 그대로 서 있었다.

아래 시뮬레이션에서 직접 조작해보자.

u_L1.00u_R0.40grid N120

play rate0.35t = 0.00 · consv 0.000 · prim 0.000

Both tracks solve the same initial jump on the same grid. Set u_R to 0.00 and the pink front stops dead while the green one keeps pace with the white dashed line. Push grid N to 320: the green error shrinks, the pink one does not — refinement never buys back a speed the scheme was never told to conserve.

위쪽 초록이 보존형, 아래쪽 분홍이 원시형이고 흰 점선이 정확한 충격파 위치다. u_R을 0.00으로 내리면 분홍 전선이 완전히 멈추고, grid N을 320까지 올려도 그 자리에서 움직이지 않는다.

검사체적에서 나온 식은 원래 flux 형태였다#

왜 flux 형태가 원본인가. 유도를 되짚으면 답이 보인다.

미소 육면체 $dx\,dy\,dz$ 를 잡고 각 면을 지나는 질량유량을 센다. 면 하나를 통과하는 양은 그 면 중심의 밀도, 면에 수직한 속도, 면적의 곱이다. 즉 $\dot m = \rho V_n A$ 다. 면 중심값은 셀 중심에서 테일러 급수로 전개하고 2차 이상을 버린다. 여섯 면을 모두 더한 뒤 $dx\,dy\,dz$ 로 나누면 연속방정식이 나온다.

\frac{\partial \rho}{\partial t} + \frac{\partial (\rho u_j)}{\partial x_j} = 0

운동량도 같은 절차다. 면을 통해 드나든 운동량, 체적력, 표면력을 더한다.

\frac{\partial (\rho u_i)}{\partial t} + \frac{\partial (\rho u_i u_j + p\,\delta_{ij})}{\partial x_j} = \frac{\partial \tau_{ij}}{\partial x_j}

여기서 $\rho$ 는 밀도, $u_i$ 는 속도 성분, $p$ 는 압력, $\tau_{ij}$ 는 점성응력 텐서(뉴턴 유체 가정에서 속도 기울기에 선형인 편차응력)다.

중요한 것은 이 식의 생김새가 아니라 출신이다. 모든 항이 "면을 통해 오간 양"으로 정의되어 있다. 발산 형태인 것은 우연이 아니라 유도 자체가 그렇게 생겼기 때문이다.

원시형은 여기서 한 단계 더 간다. 곱미분을 풀고, 연속방정식을 곱해서 빼고, $\rho$ 로 나눈다. $\rho = \text{const}$ 를 가정하면 $\partial u_i / \partial x_i = 0$ 과 함께 다음이 남는다.

\frac{\partial u_i}{\partial t} + \frac{\partial (u_i u_j)}{\partial x_j} = -\frac{1}{\rho}\frac{\partial p}{\partial x_i} + \nu\,\frac{\partial^2 u_i}{\partial x_j \partial x_j}

이 조작들은 전부 미분 가능성을 전제한다. 불연속 위에서는 전제가 없다.

나눗셈 한 번이 지워 버린 망원경 합#

이산 레벨에서 보면 더 분명하다. 보존형 유한체적법의 갱신식은 이렇다.

u_i^{n+1} = u_i^{n} - \frac{\Delta t}{\Delta x}\left(f_{i+1/2} - f_{i-1/2}\right)

셀 전체에 대해 합을 취해 보자. 내부 면 $f_{i+1/2}$ 는 셀 $i$ 에서 빠지고 셀 $i+1$ 에서 더해진다. 부호가 반대라 정확히 상쇄된다. 이것이 망원경 합(telescoping sum)이다. 남는 것은 도메인 양 끝의 flux뿐이다.

\Delta \sum_i u_i \Delta x = -\Delta t \left(f_{N+1/2} - f_{1/2}\right)

총량 변화는 경계에서 오간 양과 정확히 같다. 반올림 오차를 빼면 예외가 없다.

원시형은 여기가 깨진다. $u_i \cdot (u_i - u_{i-1})/\Delta x$ 에는 셀마다 다른 계수 $u_i$ 가 앞에 붙는다. 이웃한 두 항의 크기가 달라 상쇄되지 않는다. 남은 찌꺼기가 매 스텝 쌓인다.

코드로 잰 결과가 아래에 있다. 두 번째 리만 문제( $u_L = 1$ , $u_R = 0.4$ )에서 도메인에 실제로 들어와야 할 양은 $+0.168000$ 이다. 보존형은 소수점 여섯 자리까지 그 값을 낸다. 원시형은 $+0.152765$ , 약 9%를 잃었다.

같은 종류의 누수를 AMR 태깅 기준과 coarse-fine 리플럭싱에서 다룬 적이 있다. 그때는 격자 레벨 경계에서 면 flux가 두 개 존재하는 것이 원인이었다. 원리는 같다. 면에서 오간 양의 장부가 맞지 않으면 총량이 샌다.

Rankine–Hugoniot은 flux에만 답한다#

충격파 속도는 어디에서 오는가. 불연속을 감싸는 얇은 검사체적에 보존법칙을 적용하면 나온다.

s\,(u_L - u_R) = f(u_L) - f(u_R)

$s$ 는 불연속의 전파 속도, $f$ 는 flux다. Burgers에서 $f = u^2/2$ 이므로 다음이 된다.

s = \frac{u_L^2/2 - u_R^2/2}{u_L - u_R} = \frac{u_L + u_R}{2}

이 관계식에 들어 있는 것은 $f$ 뿐이다. $u\,\partial_x u$ 라는 표현은 등장하지 않는다. 등장할 수도 없다. 불연속에서 $\partial_x u$ 는 델타 함수이고, 거기에 점프하는 $u$ 를 곱하는 연산은 분포 이론에서 정의되지 않는다. 이것을 비보존 곱(non-conservative product)이라 부른다.

Lax–Wendroff 정리가 보증하는 것도 정확히 이 지점이다. 보존형 스킴의 수치해가 수렴하면 그 극한은 반드시 보존법칙의 약해다. 즉 Rankine–Hugoniot을 만족한다. 비보존형에는 그런 보증이 없다. Hou와 LeFloch가 보인 것은 더 나쁘다. 비보존형도 수렴은 하는데, 잘못된 속도로 수렴한다.

Python으로 잰 전파 속도와 총량#

같은 격자, 같은 CFL, 같은 초기조건으로 두 스킴을 돌렸다. 순수 표준 라이브러리만 쓴다.

def riemann_setup(nx, ul, ur, xs=0.3):
    dx = 1.0 / nx
    return dx, [ul if (i + 0.5) * dx < xs else ur for i in range(nx)]
 
def godunov_flux(a, b):
    if a > b:                                  # 충격파: 풍상 쪽을 고른다
        return 0.5 * a * a if a + b >= 0 else 0.5 * b * b
    if a >= 0:
        return 0.5 * a * a
    return 0.5 * b * b if b <= 0 else 0.0      # 음속 통과 팽창파
 
def step_conservative(u, dx, dt):              # u_t + (u^2/2)_x = 0
    n = len(u)
    f = [0.5 * u[0] ** 2] + [godunov_flux(u[i], u[i + 1]) for i in range(n - 1)] \
        + [0.5 * u[-1] ** 2]
    return [u[i] - dt / dx * (f[i + 1] - f[i]) for i in range(n)]
 
def step_primitive(u, dx, dt):                 # u_t + u u_x = 0
    n, out = len(u), []
    for i in range(n):
        im, ip = max(i - 1, 0), min(i + 1, n - 1)
        g = (u[i] - u[im]) / dx if u[i] >= 0 else (u[ip] - u[i]) / dx
        out.append(u[i] - dt * u[i] * g)
    return out
 
def shock_locate(u, dx, level):
    for i in range(1, len(u)):
        if u[i] < level <= u[i - 1]:
            return (i - 0.5) * dx + dx * (u[i - 1] - level) / (u[i - 1] - u[i])
    return float("nan")
 
def march_burgers(nx, ul, ur, tend, step):
    dx, u = riemann_setup(nx, ul, ur)
    t = 0.0
    while t < tend - 1e-12:
        dt = min(0.4 * dx / max(max(abs(v) for v in u), 1e-12), tend - t)
        u = step(u, dx, dt)
        t += dt
    return dx, u
 
T, XS = 0.4, 0.3
for ul, ur in ((1.0, 0.0), (1.0, 0.4)):
    s = 0.5 * (ul + ur)
    influx = (0.5 * ul ** 2 - 0.5 * ur ** 2) * T        # 구간에 들어와야 할 순 flux
    print("uL=%.1f uR=%.1f | Rankine-Hugoniot speed = %.3f" % (ul, ur, s))
    print("    N   conservative   primitive")
    for nx in (100, 200, 400, 800):
        v = []
        for step in (step_conservative, step_primitive):
            dx, u = march_burgers(nx, ul, ur, T, step)
            v.append((shock_locate(u, dx, s) - XS) / T)
        print("%5d      %7.4f     %7.4f" % (nx, v[0], v[1]))
    for name, step in (("conservative", step_conservative), ("primitive  ", step_primitive)):
        dx, u = march_burgers(400, ul, ur, T, step)
        dx0, u0 = riemann_setup(400, ul, ur)
        print("  N=400 %s : d(int u dx) = %+.6f  (exact %+.6f)"
              % (name, sum(u) * dx - sum(u0) * dx0, influx))
    print()

uL=1.0 uR=0.0 | Rankine-Hugoniot speed = 0.500
    N   conservative   primitive
  100       0.5006      0.0000
  200       0.5003      0.0000
  400       0.5002      0.0000
  800       0.5001      0.0000
  N=400 conservative : d(int u dx) = +0.200000  (exact +0.200000)
  N=400 primitive   : d(int u dx) = +0.000000  (exact +0.200000)
 
uL=1.0 uR=0.4 | Rankine-Hugoniot speed = 0.700
    N   conservative   primitive
  100       0.7009      0.6263
  200       0.7005      0.6330
  400       0.7002      0.6363
  800       0.7001      0.6379
  N=400 conservative : d(int u dx) = +0.168000  (exact +0.168000)
  N=400 primitive   : d(int u dx) = +0.152765  (exact +0.168000)

첫 번째 경우가 극단적이다. $u_R = 0$ 이면 불연속 오른쪽 셀에서 $u\,\partial_x u$ 가 통째로 0이다. 갱신할 것이 없으니 전선이 서지도 못하고 멈춘다. 총량 변화도 정확히 0이다. 왼쪽 경계로 들어온 $0.2$ 가 어디에도 나타나지 않는다.

격자를 조밀하게 하면 되는가#

두 번째 경우가 실무에서 더 위험하다. 원시형의 속도가 $0.6263 \to 0.6330 \to 0.6363 \to 0.6379$ 로 움직인다. 격자를 8배로 늘렸더니 값이 안정된다. 수렴하는 것처럼 보인다.

문제는 수렴 목적지다. 정답은 $0.700$ 인데 이 수열은 대략 $0.639$ 로 간다. 약 8.7% 낮다. 격자 수렴 검사(grid convergence study)를 정직하게 돌려도 이 오차는 잡히지 않는다. 세 격자에서 값이 서로 가까워지는 것을 확인하고 "수렴했다"고 적으면 끝난다.

보존형은 $0.7009 \to 0.7001$ 로 정답에 붙는다. 오차가 $\Delta x$ 에 비례해 줄어든다. 두 수열의 차이는 정확도의 차이가 아니라 푸는 방정식의 차이다.

부드러운 해만 있는 문제에서는 이 차이가 드러나지 않는다. 그래서 Taylor–Green 같은 검증 케이스만 돌린 코드는 무사히 통과한다. 불연속이 처음 생기는 순간, 그때까지 옳던 코드가 조용히 다른 물리를 풀기 시작한다. Euler 방정식의 특성곡선과 음파에서 본 특성선의 교차가 바로 그 순간이다.

그래도 원시형을 쓰는 자리 — $\rho=\text{const}$ 의 유효기간#

그렇다고 원시형이 잘못된 형태는 아니다. 비압축성 해석은 거의 전부 원시형이다. 이유가 있다.

첫째, 미지수가 줄어든다. 2차원 압축성은 $\rho, u, v, p, T$ 다섯 개를 질량·운동량 2개·에너지·상태방정식 다섯 식으로 푼다. 비압축성은 $\rho$ 를 상수로 빼고 에너지식과 상태방정식을 떼어 낸다. 남는 것은 $u, v, p$ 뿐이다.

둘째, 압력이 열역학량이 아니라 발산 구속조건의 라그랑주 승수가 된다. 그래서 압력 Poisson 방정식으로 따로 푼다. 이 구조를 Chorin의 투영법과 분수 시간전진에서 다뤘다.

셋째, 비압축성 유동에는 충격파가 없다. Rankine–Hugoniot을 지킬 불연속이 애초에 없으니 위 문제가 발생하지 않는다.

유효기간은 마하수가 정한다. 등엔트로피 관계에서 밀도는 이렇게 변한다.

\frac{\rho}{\rho_0} = \left(1 + \frac{\gamma-1}{2}M^2\right)^{-\frac{1}{\gamma-1}}

$\rho_0$ 는 정체 밀도, $\gamma$ 는 비열비, $M$ 은 마하수다. 작은 $M$ 에서 전개하면 밀도 변화가 $M^2/2$ 로 간다. $M = 0.2$ 에서 약 2%, $M = 0.3$ 에서 약 4.5%다. 흔히 쓰는 $M < 0.2$ 기준은 여기서 나온 숫자다.

exit Mach0.25drho 0.00% · speed error 0.00%

Drag exit Mach from 0.05 upward. Below 0.2 the two rows of dots stay in step and both readouts sit green — the deleted term is under 2%. Past 0.3 the pink row falls behind the green one, and the yellow dot climbs off the M²/2 dashed line: the density the incompressible model froze is now doing real work.

exit Mach를 0.05에서 끌어올리면 위 초록 점(밀도 변화 반영)과 아래 분홍 점(밀도 고정)의 간격이 벌어진다. 0.2 아래에서는 두 줄이 거의 겹치고, 0.3을 넘으면 오른쪽 노란 점이 $M^2/2$ 점선에서 떨어져 나간다.

충격파가 늦게 도착했을 때 먼저 볼 곳#

솔버가 충격파를 엉뚱한 자리에 세우면 순서대로 확인할 것이 있다.

시간전진식이 면 flux의 차분인지부터 본다. $\sum_i u_i \Delta x$ 의 변화가 경계 flux와 자릿수까지 일치해야 한다. 일치하지 않으면 다른 것을 보기 전에 이것부터다.

다음은 소스항으로 옮긴 항이다. 곡선좌표나 축대칭 항을 정리하다 보면 발산 안에 있어야 할 것이 오른쪽으로 넘어가곤 한다. 매끄러운 해에서는 아무 일도 없다가 불연속에서 속도가 틀어진다.

마지막으로 비보존 곱이 남아 있는지 본다. 다상류의 $\alpha\,\partial_x p$ 같은 항은 원리적으로 비보존이라 경로 적분 해석이 따로 필요하다. 그런 항이 있다면 격자 세분으로 해결되지 않는다는 것을 먼저 알고 있어야 한다.

세분해도 충격파 위치가 제자리인 순간, 의심할 것은 정확도가 아니라 형태다.