[논문 리뷰] β를 2로 올렸더니 계면은 14칸에 멈췄고 사인파는 1차 풍상보다 틀렸다 — Tallois 압축형 리미터

재진입체 표면에서 녹은 액막을 따라가야 했다#

대기권에 고속으로 들어오는 물체는 표면이 녹는다. 녹은 액체는 얇은 막이 되어 기체 흐름에 끌려간다. CEA-CESTA의 Tallois, Peluchon, Villedieu는 이 액막을 확산 계면(diffuse interface, 두 유체 사이를 몇 칸에 걸친 혼합 셀로 표현하는 방식) 5식 모델로 풀었다. 문제는 계면이 시간이 갈수록 번진다는 것이다. 2022년 논문은 이를 2차 MUSCL과 압축형 리미터로 막았다.

이 글은 그 리미터 한 줄, 식 (3.10)만 떼어 본다. 1D 이송 실험으로 두 가지를 확인한다. β=2는 계면을 다섯 바퀴 뒤에도 14칸에 묶어둔다. 그런데 같은 리미터가 사인파에서는 1차 풍상(upwind, 바람 불어오는 쪽 셀 값만 쓰는 1차 스킴)보다 오차가 크다. 논문이 왜 체적분율 $z$에만 β=2를 쓰는지가 이 두 숫자에 들어 있다.

번지는 계면과 리미터의 역할#

5식 모델의 체적분율은 이송 방정식을 따른다. 진짜 해에서는 0과 1 사이 계단이 모양 그대로 움직인다. 수치해는 다르다. 1차 스킴은 매 스텝 수치 확산을 더한다. 계면 폭이 $\sqrt{t}$처럼 계속 넓어진다.

2차 MUSCL(셀 안을 기울기 있는 직선으로 재구성하는 방식)은 확산을 줄인다. 대신 계단 근처에서 진동을 만든다. 그래서 기울기에 리미터 $\phi_i$를 곱한다. 논문의 재구성은 이렇다.

$\hat\alpha_{ij}$는 셀 $i$에서 이웃 $j$ 쪽 면의 재구성 값, $\kappa$는 U-MUSCL 계수(논문은 1/3), $\mathbf{x}_f$는 면 중심이다. 1D에서 $\kappa=1/3$이면 고전적인 3차 풍상 편향 재구성과 같아진다.

식 (3.10): φ를 1보다 크게 허용한다#

$|\Delta\alpha_i|_{\max}$는 모든 면 중 가장 큰 증분, $\alpha_i^{\min},\alpha_i^{\max}$는 $\alpha_i$와 이웃 기준값 $\bar\alpha_j$들의 최소·최대다. 증분이 0이면 $\phi_i=0$이다.

손잡이는 두 개다. 하나는 상한 β, 다른 하나는 기준값 $\bar\alpha_j$다.

• β=1, $\bar\alpha_j$=셀 $i$와 $j$의 거리 가중 평균 → 단조성 보존(monotonicity preserving). 면 값이 두 셀의 중간값을 넘지 못한다. minmod와 비슷하다.
• β>1, $\bar\alpha_j=\alpha_j$ → 압축형(compressive). 면 값이 이웃 셀 값까지 갈 수 있다. 기울기도 원래의 β배까지 키운다.

풍상 플럭스는 상류 셀의 하류 쪽 면 값으로 정해진다. 압축형 리미터는 그 면 값을 이웃 값 쪽으로 밀어붙인다. 계단의 앞쪽 셀은 더 많이 비우고 뒤쪽 셀은 더 빨리 채운다. 이것이 계면을 세운다.

아래 시뮬레이션에서 직접 조작해보자. 날카로운 계단이 셀 다섯 개를 지나가며, 각 셀의 재구성 직선과 오른쪽 면 값(점)이 그려진다.

계면이 가운데 셀 $i$ 위를 지날 때 주황 직선이 파란 직선보다 훨씬 가파르다. 주황 띠(이웃 값 경계)가 파란 띠(면 평균 경계)보다 넓기 때문이다. β 슬라이더를 1로 내리면 오른쪽 φ 궤적의 주황 봉우리가 1에서 잘린다.

코드: 같은 슬랩, 네 가지 리미터, 다섯 바퀴#

주기 영역 $[0,1]$에 셀 100개를 깐다. 속도는 $u=1$, CFL 0.4, 시간 적분은 논문 3.2.1절의 Heun 법이다. 초기값은 두 가지다. $0.3<x<0.7$에서 $z=1$인 슬랩과, $z=0.5+0.5\sin 2\pi x$인 매끄러운 분포다.

리미터는 네 가지를 비교한다. 1차 풍상, β=1 단조성 보존(MP), β=1에 이웃 값 경계(nbr), β=2 압축형이다. 측정값은 혼합 셀 수($0.01 < z < 0.99$), 슬랩과 사인의 L1 오차, 사인에서 $z>0.95$인 셀 수다. 정확한 사인은 14칸이 0.95를 넘는다.

import math
 
N, CFL, KAPPA = 100, 0.4, 1.0 / 3.0
h = 1.0 / N
x = [(i + 0.5) * h for i in range(N)]
 
def right_face(z, beta, midpoint_bounds):
    zf = [0.0] * N
    for i in range(N):
        zi, zm, zp = z[i], z[i - 1], z[(i + 1) % N]
        grad_h = 0.5 * (zp - zm)                       # ∇z_i · h
        # Δz_{i,j} = κ/2 (z_j - z_i) + (1-κ) ∇z_i·(x_f - x_i)   (식 3.9)
        d_r = 0.5 * KAPPA * (zp - zi) + (1 - KAPPA) * 0.5 * grad_h
        d_l = 0.5 * KAPPA * (zm - zi) - (1 - KAPPA) * 0.5 * grad_h
        if midpoint_bounds:            # β=1 단조성 보존: 면 평균값으로 경계
            bl, br = 0.5 * (zi + zm), 0.5 * (zi + zp)
        else:                          # 압축형: 이웃 셀 값으로 경계
            bl, br = zm, zp
        zmin, zmax = min(zi, bl, br), max(zi, bl, br)
        dmax = max(abs(d_l), abs(d_r))
        phi = min(beta, (zi - zmin) / dmax, (zmax - zi) / dmax) if dmax > 0 else 0.0
        zf[i] = zi + phi * d_r                         # u>0: 오른쪽 면만 필요 (식 3.10)
    return zf
 
def rhs(z, beta, mp):
    f = right_face(z, beta, mp) if beta > 0 else z     # u = 1, 풍상 플럭스
    return [-(f[i] - f[i - 1]) / h for i in range(N)]
 
def advect(z0, beta, mp, periods):
    dt = CFL * h
    z = list(z0)
    for _ in range(round(periods / dt)):
        k1 = rhs(z, beta, mp)
        z1 = [a + dt * b for a, b in zip(z, k1)]
        k2 = rhs(z1, beta, mp)
        z = [0.5 * (a + b + dt * c) for a, b, c in zip(z, z1, k2)]   # Heun
    return z
 
def l1(a, b):
    return h * sum(abs(p - q) for p, q in zip(a, b))
 
step = [1.0 if 0.3 < xi < 0.7 else 0.0 for xi in x]
wave = [0.5 + 0.5 * math.sin(2 * math.pi * xi) for xi in x]
schemes = [("1st-order", 0.0, False), ("beta=1 MP", 1.0, True),
           ("beta=1 BJ", 1.0, False), ("beta=2", 2.0, False)]
 
print("scheme      periods  mixture  L1(step)  L1(sine)  sine_flat")
for name, beta, mp in schemes:
    for periods in (1, 5):
        zs = advect(step, beta, mp, periods)
        zw = advect(wave, beta, mp, periods)
        mix = sum(1 for v in zs if 0.01 < v < 0.99)
        flat = sum(1 for v in zw if v > 0.95)          # 꼭대기가 평평해진 셀 수
        print(f"{name:10s} {periods:6d} {mix:8d} {l1(zs, step):9.4f}"
              f" {l1(zw, wave):9.4f} {flat:9d}")
print("exact sine: cells with z>0.95 =", sum(1 for v in wave if v > 0.95))

scheme      periods  mixture  L1(step)  L1(sine)  sine_flat
1st-order       1       87    0.1595    0.0570         0
1st-order       5      100    0.3427    0.1997         0
beta=1 MP       1       48    0.0687    0.0054        14
beta=1 MP       5       81    0.1242    0.0228         0
beta=1 BJ       1       18    0.0337    0.0011        14
beta=1 BJ       5       26    0.0496    0.0034        15
beta=2          1       14    0.0170    0.0616        31
beta=2          5       14    0.0170    0.1636        46
exact sine: cells with z>0.95 = 14

위 코드를 브라우저에서 직접 돌려보자. 네 리미터가 슬랩과 사인을 동시에 다섯 바퀴 이송하고, 오른쪽에 혼합 셀 수와 사인 L1 오차의 이력이 쌓인다.

기본값에서 1바퀴 시점의 숫자가 출력 블록과 같다. 주황(β)을 1.5로 내리면 사인 오차가 줄고 혼합 셀은 늘어난다. 셀 수 N을 200으로 올리면 1차 풍상의 혼합 셀은 늘지만 β=2는 여전히 몇 칸에 머문다.

14칸과 0.0616이 함께 말하는 것#

슬랩 쪽 숫자부터 보자. 1차 풍상은 1바퀴 만에 87칸이 혼합 셀이 됐다. 5바퀴 뒤에는 100칸 전부다. 계단이 사라졌다. β=1 MP도 48칸에서 81칸으로 계속 번진다. β=2는 1바퀴와 5바퀴 모두 14칸이다. L1 오차도 0.0170으로 같다.

14칸이 멈춘 것은 우연이 아니다. 압축형 리미터는 계면을 일정한 모양의 진행파로 만든다. 수치 확산과 압축이 균형을 이룬 폭이다. 계면 하나당 약 7칸이고, 0.05~0.95 구간은 3칸 남짓이다.

이제 사인 쪽이다. β=2의 1바퀴 L1 오차는 0.0616이다. 1차 풍상의 0.0570보다 크다. 2차 정확도 스킴이 1차보다 틀린 셈이다. 마지막 열이 이유를 보여준다. 0.95를 넘는 셀이 14칸에서 31칸, 5바퀴 뒤 46칸으로 늘었다. 봉우리가 평평해지고 비탈이 가팔라졌다. 사인이 사각파로 변한 것이다. superbee가 사인을 깎는 것과 같은 현상이다.

β=1 nbr도 눈여겨볼 만하다. 슬랩 혼합 셀은 18→26칸으로 천천히 번진다. 사인 오차는 0.0011로 넷 중 가장 작다. 경계를 이웃 값으로 넓히는 것만으로도 MP보다 훨씬 덜 번진다. 대신 MP가 주는 면 값 순서 보장(식 3.11)은 없다.

그래서 논문은 z에만 β=2를 준다#

논문 3.1.3절의 배치는 이 결과와 정확히 맞는다.

• 음향 단계: $(u,p)$를 β=1로 재구성한다. 이웃은 면 이웃 5점이다.
• 이송 단계: $(\rho_1,\rho_2,u,p)$는 β=1, 체적분율 $z$만 β=2다. 이웃은 꼭짓점 이웃까지 9점이다.

밀도·속도·압력은 매끄러운 부분이 많다. 거기에 β=2를 걸면 위의 사인처럼 계단이 생긴다. 게다가 음향 단계는 암시적이다. Toro–Vázquez 분할 글에서 본 것처럼 강성이 큰 압력 부분계에서는 견고함이 우선이다. 반면 $z$는 원래 계단이다. 사각파로 깎여도 잃을 것이 없다.

재구성 변수가 원시 변수라는 점도 중요하다. 보존 변수는 식 (3.12)로 다시 만든다. $z$를 날카롭게 하면서 압력을 건드리지 않아야 계면의 압력 진동을 피한다. 이 5식 모델의 $z$ 식에 붙은 $K\nabla\cdot u$ 항은 Murrone–Guillard 글에서 다뤘다.

Zalesak 원판과 1D 실험의 거리#

논문은 Zalesak 원판(홈이 파인 원판을 한 바퀴 회전시키는 시험)으로 리미터를 검증한다. 256×256 격자, CFL 0.4다. $\kappa=1/3$이 $\kappa=0$보다 홈과 모서리를 잘 살렸다고 보고한다.

1D 실험이 다루지 못한 것도 있다. 비정렬 격자에서는 $|\Delta\alpha_i|_{\max}$가 여러 면의 최대값이다. 한 면 때문에 다른 면의 기울기까지 잘린다. 또 2D 회전 이송에서는 압축형 리미터가 격자 방향으로 계면을 계단처럼 꺾는 경향이 있다. 원판의 둥근 가장자리가 그 시험대다.

압축형 리미터를 다른 변수에 옮기기 전에#

압축형 리미터는 "계단이어야 하는 변수"에만 쓴다. 내 코드에 β>1을 넣기 전에 매끄러운 사인 하나를 한 바퀴 돌려보자. L1 오차가 1차 풍상보다 크다면 그 변수에는 β=1이 맞다. 혼합 셀 수가 바퀴마다 늘어난다면 $z$의 β가 부족하다는 신호다.