[논문 리뷰] 둘 중 더 작은 진폭을 골라라 — Deng(2018) MUSCL-THINC-BVD 인터페이스 재구성

한 칸 순수 점프는 차분 곱이 0이라, THINC는 경계 변동 비교 전에 탈락한다.

도쿄공업대 Xiao 그룹이 2018년에 두께를 재 봤다. 압축성 이상 유동, 1000스텝. 한 칸이던 점프가 여덟 칸이 됐다. WENO도 못 막았다. 스텝마다 소산이 작아 보여도, 시간이 쌓이면 계면은 무너진다.

그들이 들고 온 재구성이 MUSCL-THINC-BVD다. 셀마다 후보를 둘 만든다. 경계 점프가 더 작은 쪽만 남긴다.

β 슬라이더를 1.4에서 2.2로 올리면 점프 칸의 tanh가 좁아지고, 가우시안 어깨는 톱니로 반응한다.

N = 80β = 1.60CFL = 0.40

MUSCL THINC BVDt = 0.00THINC picked: 0% of cells

사각파 가장자리에만 청록 막대가 켜지는지 본다. 1.6 근처에서 그 막대가 더 늘어나지 않는지도 같이 본다.

한 페이지 요약#

저자 / 저널: Deng, Inaba, Xie, Shyue, Xiao. Journal of Computational Physics 371 (2018) 945–966.
문제: 5-방정식 모델로 푸는 압축성 이상 유동에서 물질 계면이 시간에 따라 번진다.
제안: 매 셀에서 MUSCL(매끈한 구간)과 THINC(점프)를 같이 만들고, 경계 점프가 작은 쪽을 남기는 BVD(Boundary Variation Diminishing).
다른 점: 후처리 anti-diffusion이나 인공 압축이 없다. 부피분율과 나머지 보존 변수에 같은 규칙을 걸어, 변수 사이 정합이 따라온다.

충돌하는 두 요구#

압축성 이상 유동은 두 가지를 동시에 원한다. 매끈한 곳은 정확하고 소산이 적어야 한다. 점프, 그러니까 계면과 충격에서는 단조성과 얇은 두께가 필요하다. 함수 하나가 둘 다 잘하기는 어렵다.

MUSCL은 구간마다 직선이다. 단조성은 지키지만 기울기는 1차까지다. 계면을 지날 때마다 조금 깎인다. 1000스텝이면 한 칸이 여덟 칸이 되는 누적이다. THINC는 셀 안에 tanh를 맞춰 계단을 한두 칸에 가둔다. 매끈한 곳에 그대로 쓰면 가짜 계단이 생긴다.

논문의 출발은 결정이다. 둘을 한 함수로 섞지 않는다. 셀마다 하나를 고른다.

MUSCL — 신뢰할 수 있는 디퓨저#

기준 후보는 minmod 기울기 제한을 쓴 MUSCL이다.

\tilde{q}_i^{\,\text{MUSCL}}(x) = \bar{q}_i + \sigma_i (x - x_i), \quad \sigma_i = \text{minmod}(\bar{q}_i - \bar{q}_{i-1},\; \bar{q}_{i+1} - \bar{q}_i)

$\bar q_i$ 는 셀 평균이다. $\sigma_i$ 는 중심에서의 기울기다. minmod는 두 차분의 부호가 같으면 절대값이 작은 쪽을, 다르면 0을 준다.

경계값은 이렇게 읽는다.

q^{L,\text{MUSCL}}_{i+1/2} = \bar{q}_i + \tfrac{1}{2}\sigma_i, \quad q^{R,\text{MUSCL}}_{i-1/2} = \bar{q}_i - \tfrac{1}{2}\sigma_i

진동은 없다. 점프 구간에서 두께만 천천히 벌어진다.

THINC — 점프를 흉내내는 단조함수#

THINC(Tangent of Hyperbola for INterface Capturing)는 셀 안의 점프를 쌍곡탄젠트로 맞춘다.

\tilde{q}_i^{\,\text{THINC}}(x) = \bar{q}_{\min} + \tfrac{\bar{q}_{\max}}{2}\left[1 + \theta \tanh\!\left(\beta\!\left(\tfrac{x - x_{i-1/2}}{\Delta x} - \tilde{x}_i\right)\right)\right]

$\bar q_{\min}$ 과 $\bar q_{\max}$ 는 이웃 평균에서 잡은 최소와 진폭이다. $\theta$ 는 점프 방향이다. $\beta$ 는 두께다. $\tilde x_i$ 는 셀 평균이 유지되도록 푸는 점프 중심이다.

$\beta$ 는 1.4에서 2.0 사이가 안정적이다. 표준은 1.6이다. 더 크면 점프가 한 칸에 갇혀 보기엔 좋다. 매끈한 곳에 들어가면 톱니가 된다. tanh를 충격 쪽에 붙인 기록은 TENO와 THINC를 한 스텐실에 겹친 재구성에 따로 있다. 거기는 후보를 고르는 대신 가중으로 섞는다.

BVD — 두 후보 사이의 변동 재기#

후보가 둘이면, 이웃 재구성과 경계에서 만날 때 어느 쪽 점프가 더 작은지 재면 된다. 셀 $i$ 의 총 경계 변동(Total Boundary Variation, TBV)은 이렇게 정의된다.

\text{TBV}_i^{P} = |q^{L,\text{MUSCL}}_{i-1/2} - q^{R,P}_{i-1/2}| + |q^{L,P}_{i+1/2} - q^{R,\text{MUSCL}}_{i+1/2}|

$P$ 는 MUSCL 또는 THINC다. 양쪽 이웃은 MUSCL로 고정한다. 가운데 셀만 바꾼다.

선택 규칙은 이렇다.

\tilde{q}_i^{\,\text{BVD}}(x) = \begin{cases} \tilde{q}_i^{\,\text{THINC}}(x) & \text{if } \delta < C_i < 1 - \delta,\ (\bar{q}_{i+1}-\bar{q}_i)(\bar{q}_i-\bar{q}_{i-1}) > 0,\ \text{TBV}^{\text{THINC}}_i < \text{TBV}^{\text{MUSCL}}_i \\ \tilde{q}_i^{\,\text{MUSCL}}(x) & \text{otherwise} \end{cases}

$C_i$ 는 점프가 셀의 어디쯤 있는지다. $\delta$ 는 약 $10^{-4}$ 다.

선택식의 첫 관문은 차분 곱이다. 이웃 차분 둘을 곱해 0보다 커야 한다. 한쪽이 0이면 곱도 0이라 관문을 넘지 못한다. 0과 1이 한 칸에 붙은 계단이 그 모양이다. 점프 직전 셀은 오른쪽 차분만 있고 왼쪽 차분은 없다. 총 경계 변동을 재기 전에 MUSCL로 돌아간다.

두 번째 관문은 위치 척도다. $\delta$ 는 $10^{-4}$ 근처다. 이웃이 잡은 진폭의 그 비율 안으로 셀 평균이 붙으면 THINC 자격이 없다. 평탄한 마루를 잘라 내는 장치다.

5-방정식 모델 위의 일관성#

재구성 비교에서 멈추지 않는 지점이 여기다. 5-방정식 모델은 부피분율 $\alpha_1$ , 상 밀도 $\alpha_1\rho_1$ , $\alpha_2\rho_2$ , 운동량 $\rho u$ , 총에너지 $E$ 를 같이 푼다. 다섯 변수에 같은 BVD 규칙을 건다.

예전 방법은 부피분율만 날카롭게 잡고, 나머지는 압력 진동이 없도록 따로 고쳤다. BVD는 다섯 변수가 같은 칸에서 같은 결정을 따르게 한다. 점프 칸이면 다섯 모두 THINC다. 매끈하면 다섯 모두 MUSCL이다. 변수 사이 정합이 따라오고, anti-diffusion 단계를 따로 두지 않는다.

NumPy로 본 BVD 결정 트리#

import numpy as np
 
def muscl_minmod_edges(q):
    """ q[i-1], q[i], q[i+1] 을 받아 MUSCL 좌/우 경계값을 돌려준다 """
    qm, q0, qp = q
    a, b = q0 - qm, qp - q0
    if a * b <= 0.0:
        slope = 0.0
    else:
        slope = np.sign(a) * min(abs(a), abs(b))
    return q0 - 0.5 * slope, q0 + 0.5 * slope  # qL, qR
 
def thinc_jump_edges(q, beta=1.6, eps=1e-20):
    """ tanh 점프를 fit한 좌/우 경계값 """
    qm, q0, qp = q
    qmin = min(qm, qp)
    qmax = max(qm, qp) - qmin
    if qmax < 1e-12:
        return q0, q0
    theta = np.sign(qp - qm)
    C = (q0 - qmin + eps) / (qmax + eps)
    if C <= 1e-6 or C >= 1 - 1e-6:
        return q0, q0
    B = np.exp(theta * beta * (2.0 * C - 1.0))
    A = (B / np.cosh(beta) - 1.0) / np.tanh(beta)
    qR = qmin + 0.5 * qmax * (1.0 + theta * A)
    num = 1.0 + theta * A * np.tanh(beta) + theta * np.tanh(beta)
    den = 1.0 + A * np.tanh(beta)
    qL = qmin + 0.5 * qmax * (num / den)
    return qL, qR
 
def bvd_decide(q_window, beta=1.6, delta=1e-4):
    """
    q_window: 길이 5 배열 [q[i-2], q[i-1], q[i], q[i+1], q[i+2]]
    반환: (qL_i, qR_i, picked_thinc)
    """
    qm2, qm1, q0, qp1, qp2 = q_window
    monotone = (qp1 - q0) * (q0 - qm1) > 0.0
    qmin = min(qm1, qp1)
    qmax = max(qm1, qp1) - qmin
    C = 0.5 if qmax < 1e-12 else (q0 - qmin) / qmax
    eligible = monotone and (delta < C < 1.0 - delta)
 
    qL_M, qR_M = muscl_minmod_edges([qm1, q0, qp1])
    if not eligible:
        return qL_M, qR_M, False
 
    qL_T, qR_T = thinc_jump_edges([qm1, q0, qp1], beta=beta)
    # 양쪽 이웃은 MUSCL로 고정한 TBV 비교
    _, qR_left = muscl_minmod_edges([qm2, qm1, q0])
    qL_right, _ = muscl_minmod_edges([q0, qp1, qp2])
    tbv_M = abs(qR_left - qL_M) + abs(qR_M - qL_right)
    tbv_T = abs(qR_left - qL_T) + abs(qR_T - qL_right)
    if tbv_T < tbv_M:
        return qL_T, qR_T, True
    return qL_M, qR_M, False
 
# 시험: 매끈 + 점프 혼합 프로파일
N = 80
x = (np.arange(N) + 0.5) / N
q = np.where((x >= 0.6) & (x <= 0.8), 1.0, 0.0)
q += 0.5 * np.exp(-((x - 0.25) / 0.04) ** 2)
 
picks = 0
for i in range(2, N - 2):
    _, _, used_thinc = bvd_decide(q[i-2:i+3])
    picks += int(used_thinc)
print(f"THINC 채택률: {picks}/{N-4}")
# 출력 예: THINC 채택률: 4/76  (점프 양쪽 두 셀씩만)

스니펫 끝 주석은 채택 4/76을 적는다. 바로 위 관문과 그 숫자는 같은 모양을 가리키지 않는다. 코드의 사각파는 조건문으로 0 또는 1이다. 한 칸 점프라 차분 하나가 0이다. 주석의 네 칸은 기울기가 셀 두세 개에 걸쳐 같은 부호로 올라가는 램프를 세야 나온다. 1000스텝 뒤의 여덟 칸은 이미 그 램프다. 두께가 한 칸으로 남아 있는 동안 이 분기는 THINC를 부르지 않는다. 번진 뒤에야 부른다.

직선에서는 반대다. 이웃 차분이 같으면 minmod 기울기는 그 차분이다. 왼쪽 면값은 왼쪽 이웃의 오른쪽 면값과 만나고, 오른쪽도 같다. MUSCL 쪽 총 경계 변동은 0이다. THINC의 tanh는 직선에 굽힘을 넣는다. 면에 틈이 생기고 변동이 커진다. 부등호는 MUSCL 쪽이다. 가우시안이 셀보다 넓으면 국소적으로 이 직선에 가깝다. 매끈한 구간이 후보 교체를 거부하는 이유다.

매끈한 가우시안은 국소적으로 직선에 가깝다. 그래서 교체가 거절된다. 사각파가 한 칸으로 붙어 있는 동안도 관문은 열리지 않는다.

셀마다 누가 이기는가 (직접 보기)#

앞의 이류는 세 재구성을 같이 밀고, BVD가 THINC를 고른 칸을 막대로 칠한다. $\beta$ 를 1.4에서 2.2로 올리면 그 칸의 점프만 좁아진다. 너무 크면 가우시안 어깨에 톱니가 보인다. 1.6에서 막대 수가 더 늘지 않는 자리가 표준값이다. 사각파가 조금 번진 뒤에야 가장자리에 청록이 붙는다. 두께가 한 칸일 때는 막대가 비어 있다.

어디까지 믿을 수 있을까#

수렴은 2차에 못 미친다. 후보의 한쪽이 MUSCL이라, 매끈한 구간의 정확도는 MUSCL에서 멈춘다. WENO 같은 고차 후보로 바꾸면 5차도 가능하지만, 이 논문은 단순한 쪽을 택했다.

TBV는 이웃이 MUSCL이라고 가정한다. 양쪽이 같이 THINC로 가는 경우는 한 번의 쓸기로 무시한다. 논문 4절의 그림은 그 단순화로도 괜찮았다. 더 일반적인 형태는 뒤로 미뤄졌다.

다차원으로 가면 $\beta$ 가 면 법선 방향으로 갈린다. 계면 법선은 Young 알고리즘으로 추정한다. 꼭짓점과 꺾임에서는 법선이 틀어져 점프가 한쪽으로 기운다.

시간 적분도 따로다. 논문은 SSPRK3를 쓴다. BVD 판정은 Runge–Kutta 단계마다 다시 해야 한다. 단계 사이에서 판정이 흔들리면 계면이 떨린다. 그 떨림은 후속 연구가 다시 집었다.

한 칸 점프에서 고르기가 멈추는 자리#

한 칸으로 붙어 있는 동안 이 분기는 THINC를 부르지 않는다. 두께가 두세 칸으로 번진 뒤에야 차분 곱이 양수가 된다. 계면을 한 칸에 붙잡는 일과, 붙잡은 뒤에야 고를 수 있는 일이 서로 어긋난다.

그 어긋남이 시간 간격 쪽으로 새는지가 다음 문제다. 계면을 얇게 유지해도 이류의 상한은 남는다. 모세관 항을 걷어낸 뒤에도 CFL 0.05가 남은 계면 이류는 그 상한이 어디 있는지 묻는다. BVD가 점프를 한 칸에 가둔다면, 그 칸을 옮기는 시간 간격은 어디까지 커질 수 있는가.