방향을 쪼개지 않고 이류를 푼다 — Bell-Colella-Glaz 비분할 업윈드

1980년대 후반, 로런스 버클리 연구소의 John Bell과 Phillip Colella는 폭발 시뮬레이션에서 이상한 자국을 발견했다. 격자 축에 비스듬히 놓인 화염면이 시간이 갈수록 사선 방향으로 번지고 비틀렸다. 범인은 물리가 아니라 이산화였다. 한 방향씩 번갈아 푸는 방향분할(dimensional splitting)이 대각선 이류를 일그러뜨린 것이다. 이 글은 그 해법으로 나온 Bell-Colella-Glaz(BCG) 비분할 2차 업윈드 스킴을 수식부터 numpy 구현까지 따라간다. 읽고 나면, 왜 회전하는 스칼라를 풀 때 방향을 쪼개면 안 되는지, 그리고 횡방향 미분 한 항이 어떻게 그 문제를 메우는지 손으로 만져보게 된다.

다루는 방정식은 가장 단순한 보존형 이류식이다. 발산이 0인(divergence-free) 속도장 $\mathbf{u}$ 가 스칼라 $s$ 를 실어 나른다.

\frac{\partial s}{\partial t} + \nabla\cdot(\mathbf{u}\,s) = 0

여기서 $s$ 는 농도·온도 같은 수동 스칼라, $\mathbf{u}$ 는 미리 주어진 속도장이다. 유한체적(FV)으로 셀 평균을 갱신하면 다음 형태가 된다.

s_{ij}^{n+1} = s_{ij}^{n} - \frac{\Delta t}{\Delta x}\left(F^x_{i+\frac12,j} - F^x_{i-\frac12,j}\right) - \frac{\Delta t}{\Delta y}\left(F^y_{i,j+\frac12} - F^y_{i,j-\frac12}\right)

$F^x = u\,s^{n+1/2}$ 는 face에서의 시간중심(time-centered) 플럭스다. 모든 정확도는 결국 face 위 값 $s^{n+1/2}$ 를 얼마나 잘 추정하느냐에 달려 있다.

한 방향씩 쪼개면 회전이 비틀린다#

방향분할은 2D 문제를 x 방향 1D 업데이트 다음 y 방향 1D 업데이트로 쪼갠다. 구현은 쉽다. 1D 솔버 하나면 끝이다. 문제는 대각선 흐름이다.

스칼라가 $45^\circ$ 로 흐르면, x를 먼저 풀고 y를 푸는 순서와 그 반대 순서가 다른 답을 낸다. 두 번의 1D 스텝 사이에 끼는 횡방향 결합(cross term)을 빠뜨리기 때문이다. 회전 유동에서는 이 오차가 매 스텝 누적되어 원판이 마름모로 일그러진다. Strang 분할로 순서를 번갈아도 2차 오차의 방향 비대칭만 줄어들 뿐, 근본 원인은 남는다.

BCG의 처방은 간단하다. x도 y도 동시에 푼다. 한 스텝 안에서 모든 face 값을 같은 시간 $t^{n+1/2}$ 로 예측하고, 그 예측에 횡방향 변화를 미리 반영한다.

시간을 공간으로 바꾸는 예측자#

핵심 아이디어는 face 값을 Taylor 전개로 시간·공간 양쪽으로 반 스텝 밀어내는 것이다. 셀 $i$ 의 오른쪽 face에서 좌측 상태를 보면,

s^{n+1/2}_{i+\frac12,\,L} = s^n_{ij} + \frac{\Delta x}{2}\,\partial_x s + \frac{\Delta t}{2}\,\partial_t s

여기서 시간 미분을 그대로 두면 쓸모가 없다. 이류 방정식 $\partial_t s = -u\,\partial_x s$ 를 대입해 시간 미분을 공간 미분으로 바꾼다(Cauchy-Kowalewski 트릭).

s^{n+1/2}_{i+\frac12,\,L} = s^n_{ij} + \frac{1}{2}\left(1 - \frac{\Delta t}{\Delta x}\,u_{i+\frac12,j}\right)\Delta_x s_{ij}

$\Delta_x s_{ij}$ 는 셀 $i$ 의 제한된 기울기(limited slope), $u_{i+1/2,j}$ 는 face의 이류 속도다. 괄호 안 $(1 - \mathrm{CFL})$ 항이 시간 전진을 담는다. CFL이 1에 가까우면 좌측 상태의 기여가 줄고, 0이면 단순한 공간 외삽이 된다.

오른쪽 셀에서 본 우측 상태도 대칭으로 만든다.

s^{n+1/2}_{i+\frac12,\,R} = s^n_{i+1,j} - \frac{1}{2}\left(1 + \frac{\Delta t}{\Delta x}\,u_{i+\frac12,j}\right)\Delta_x s_{i+1,j}

두 후보 중 하나를 고르는 것은 Godunov의 업윈드 규칙이다. face 속도의 부호가 바람을 가른다.

s^{n+1/2}_{i+\frac12} = \begin{cases} s_{i+\frac12,\,L}, & u_{i+\frac12} > 0 \\[2pt] s_{i+\frac12,\,R}, & u_{i+\frac12} < 0 \end{cases}

아래 시뮬레이션에서 직접 조작해보자. 1차원 셀 평균(계단)에서 제한된 기울기와, 그 위에서 예측된 좌측 edge 상태(주황 점)가 어떻게 정해지는지 보여준다.

CFL = 0.50

Cyan = limited PLM slope inside each cell · Amber dot = predicted left edge state at n+1/2 (advecting velocity u>0). Pick none near the jump to see the slope overshoot the neighbor values.

CFL을 0.9까지 올리면 주황 점이 셀 평균 쪽으로 끌려온다. 시간 전진이 그만큼 크기 때문이다. 점프 근처에서 limiter를 none으로 두면 PLM 직선이 이웃 값을 넘어 솟구치는 것을 볼 수 있다.

횡방향 미분이 모서리를 잇는다#

방향분할이 빠뜨린 결합을 BCG는 예측자 안에 직접 집어넣는다. x-face의 좌측 상태에 횡방향 플럭스 발산을 반 스텝만큼 빼준다.

s^{n+1/2}_{i+\frac12,\,L} \;\mathrel{-}=\; \frac{\Delta t}{2}\,\frac{\partial (v\,s)}{\partial y}\bigg|_{ij}

$v$ 는 y 방향 속도, $\partial_y(v s)$ 는 셀 $ij$ 에서의 횡방향 이류 변화다. 이 항을 추정하려면 먼저 y-face 위의 1D 예측 상태 $\hat s^y$ 를 구해 두어야 한다. 즉 BCG는 두 단계다. 먼저 각 방향의 1D 업윈드 상태를 예측하고, 그 예측을 횡방향 보정에 재사용한 뒤, 같은 업윈드 규칙으로 최종 face 값을 고른다. 이 "모서리 결합(corner coupling)"이 대각선 이류를 비틀림 없이 옮긴다.

기울기에 재갈을 물리다 — minmod와 superbee#

예측자의 정확도는 기울기 $\Delta_x s$ 에서 나온다. 중심차분 기울기를 그냥 쓰면 2차 정확도지만, 불연속 근처에서 진동(overshoot)을 만든다. 그래서 경사 제한자(slope limiter)로 재갈을 물린다. 가장 보수적인 minmod는 두 한쪽 차분 중 작은 쪽을, 부호가 다르면 0을 택한다.

\Delta_x s_{ij} = \operatorname{minmod}\!\left(s_{ij}-s_{i-1,j},\; s_{i+1,j}-s_{ij}\right)

superbee는 같은 단조성을 지키면서 기울기를 최대한 살려 접촉면을 더 날카롭게 끈다. 위 1D 데모에서 minmod와 superbee를 번갈아 보면, superbee 쪽 PLM 직선이 더 가파른 것을 확인할 수 있다.

numpy로 돌리는 회전 원판#

이제 2D 비분할 BCG를 돌려본다. 토이 문제는 강체 회전(solid-body rotation)이다. 중심을 도는 발산-0 속도장에 스칼라 원판을 띄우고, 1차 업윈드와 BCG를 같은 격자에서 비교한다.

import numpy as np
 
N = 64
dx = 1.0 / N
omega = 2 * np.pi                      # 단위 시간당 1회전
 
def face_velocities():
    j = np.arange(N)
    i = np.arange(N)
    u = -omega * ((j + 0.5) * dx - 0.5)        # x-face 속도 (y에만 의존)
    u = np.broadcast_to(u[None, :], (N, N)).copy()
    v = omega * ((i + 0.5) * dx - 0.5)         # y-face 속도 (x에만 의존)
    v = np.broadcast_to(v[:, None], (N, N)).copy()
    return u, v
 
def minmod_pair(a, b):                          # 단조성 보존 기울기
    return np.where(a * b > 0, np.where(np.abs(a) < np.abs(b), a, b), 0.0)
 
def limited_slopes(s):
    slx = minmod_pair(s - np.roll(s, 1, 0), np.roll(s, -1, 0) - s)
    sly = minmod_pair(s - np.roll(s, 1, 1), np.roll(s, -1, 1) - s)
    return slx, sly
 
def donor_fluxes(s, u, v):                       # 1차 업윈드 (도너 셀)
    Fx = np.where(u > 0, u * np.roll(s, 1, 0), u * s)
    Fy = np.where(v > 0, v * np.roll(s, 1, 1), v * s)
    return Fx, Fy
 
def bcg_fluxes(s, u, v, dt):                     # 비분할 2차 예측자
    slx, sly = limited_slopes(s)
    # 1) 법선 방향 1D 예측 상태
    sLx = np.roll(s, 1, 0) + 0.5 * (1 - dt/dx * u) * np.roll(slx, 1, 0)
    sRx = s - 0.5 * (1 + dt/dx * u) * slx
    hatX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sBy = np.roll(s, 1, 1) + 0.5 * (1 - dt/dx * v) * np.roll(sly, 1, 1)
    sTy = s - 0.5 * (1 + dt/dx * v) * sly
    hatY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    # 2) 횡방향(모서리) 보정 — 예측 상태를 재사용
    divVy = (np.roll(v * hatY, -1, 1) - v * hatY) / dx
    divUx = (np.roll(u * hatX, -1, 0) - u * hatX) / dx
    sLx -= 0.5 * dt * np.roll(divVy, 1, 0)
    sRx -= 0.5 * dt * divVy
    sBy -= 0.5 * dt * np.roll(divUx, 1, 1)
    sTy -= 0.5 * dt * divUx
    # 3) 같은 업윈드 규칙으로 최종 face 값 선택
    sX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    return u * sX, v * sY
 
def advance(s, u, v, dt, scheme):
    Fx, Fy = bcg_fluxes(s, u, v, dt) if scheme == 'bcg' else donor_fluxes(s, u, v)
    return s - dt * ((np.roll(Fx, -1, 0) - Fx) / dx + (np.roll(Fy, -1, 1) - Fy) / dx)
 
def make_disk():
    X, Y = np.meshgrid((np.arange(N)+0.5)*dx, (np.arange(N)+0.5)*dx, indexing='ij')
    return ((X - 0.5)**2 + (Y - 0.78)**2 < 0.13**2).astype(float)
 
u, v = face_velocities()
Umax = omega * 0.5 * np.sqrt(2)
dt = 0.6 * dx / Umax                    # CFL = 0.6
steps = int(round(1.0 / dt))            # 정확히 한 바퀴
 
for scheme in ('upwind1', 'bcg'):
    s = make_disk()
    mass0 = s.sum()
    for _ in range(steps):
        s = advance(s, u, v, dt, scheme)
    print(f"{scheme:8s}  peak={s.max():.3f}  "
          f"mass_err={abs(s.sum()-mass0)/mass0:.2e}")

출력은 다음과 비슷하다.

upwind1   peak=0.349  mass_err=3.1e-16
bcg       peak=0.806  mass_err=2.8e-16

두 스킴 모두 질량은 기계 정밀도까지 보존한다. 보존형 플럭스 차분이라 당연하다. 갈리는 것은 봉우리(peak)다. 1차 업윈드는 한 바퀴 만에 원판을 절반 이하로 뭉개고, BCG는 원래 높이의 80%를 지킨다.

두 바퀴를 돈 뒤 남은 것#

수치만으로는 감이 안 온다. 아래 캔버스에서 직접 회전시켜보자. 스킴 버튼으로 1차 업윈드와 BCG를 전환하고, CFL을 바꿔가며 원판이 어떻게 변형되는지 관찰한다.

CFL = 0.60

A scalar disk rotates rigidly about the center. Switch to 1st-order upwind and watch the disk melt into a ring within one revolution; BCG unsplit keeps its edge far longer.

1st-order upwind로 두면 원판이 한 바퀴도 못 돌아 고리처럼 퍼지고 가장자리가 흐려진다. BCG unsplit으로 바꾸면 경계가 훨씬 오래 살아남는다. CFL을 0.9로 올리면 BCG도 약간 거칠어지지만, 1차 업윈드와의 격차는 그대로다. 수치 확산(numerical diffusion)이 격자 해상도가 아니라 스킴 차수에서 갈린다는 사실이 눈에 들어온다.

옮길 때 새는 곳#

CFL은 분할 스킴보다 빡빡하다. 비분할 2D는 대각선 정보 전파 때문에 안정 한계가 1차원보다 작다. 실무에서는 $\mathrm{CFL}\lesssim 0.8$ 로 잡고 시작한다.
속도장은 반드시 face에 둔다. 셀 중심 속도를 face로 보간하면 발산-0이 깨지고, 보존된 줄 알았던 질량이 샌다. MAC 배치(staggered)가 정답이다.
횡방향 항을 빼먹지 말 것. 법선 예측만 구현하면 그냥 차원분할 MUSCL이다. 코드가 돌긴 하지만 회전 비틀림은 그대로 남는다.
limiter는 극값에서 0이 된다. 봉우리 꼭대기에서 기울기가 0으로 깎이면서 국소 1차로 떨어진다. 이게 봉우리가 조금씩 낮아지는 이유다. 더 살리려면 극값 보존 limiter를 검토한다.

한 줄로 남기는 것#

BCG의 정수는 "예측-보정"이다. 시간 미분을 이류로 바꿔 face 값을 반 스텝 앞으로 예측하고, 그 예측을 횡방향 보정에 재사용한 뒤, Godunov 업윈드로 최종 값을 고른다. 방향을 쪼개지 않기에 회전이 비틀리지 않고, limiter 덕에 불연속에서 진동하지 않는다. 회전·전단이 섞인 유동에서 스칼라를 옮길 일이 있다면, 1차 업윈드의 유혹을 한 번 더 의심해보자.