번진 계면을 매 스텝 되돌리다 — 안티확산 계면 샤프닝(IST) 직접 구현

얇게만 만들면 칸 수가 다시 흔들린다. 후처리로 ε/a를 고정해야 충격 옆에서도 같은 폭이 남는다.

shock-bubble 문제를 5-방정식 모델로 돌렸다. 200스텝 뒤, helium 기포의 경계는 6칸에서 20칸으로 번져 있었다. 물리량은 맞는데 계면이 죽처럼 뭉개진다. 이게 확산 계면법(diffuse interface method·계면을 유한한 폭으로 번지게 두는 접근)의 숙명이다. Nguyen 등(2021)은 이 번짐을 매 시간스텝 뒤 후처리 한 번으로 되돌린다. 오늘은 그 안티확산 기법을 직접 짜서, 왜 계면 두께가 상수로 고정되는지 눈으로 확인한다.

numerical diffusion을 켜면 6칸 경계가 칸을 삼키고, sharpening (IST)을 켜면 그 퍼짐이 ε/a 근처에서 멈춘다.

compression a: 1.00regularization ε: 0.0040

measured band: 0 cells

target ≈ ε/a: 2 cells

steady tanh thickness scales like ε/a

노란 띠가 목표 두께다. compression a를 올리면 띠가 줄고, regularization ε를 올리면 다시 두꺼워진다. 측정 칸 수가 그 비를 따라가는지만 보면 된다.

논문 정보#

제목: Numerical modeling of multiphase compressible flows with the presence of shock waves using an interface-sharpening five-equation model
저자: Van-Tu Nguyen, Thanh-Hoang Phan, Warn-Gyu Park (Pusan National University)
출처: International Journal of Multiphase Flow, 2021, 103542
DOI: 10.1016/j.ijmultiphaseflow.2020.103542
한 줄 요약: 5-방정식 압축성 2상 모델에 계면 샤프닝(IST)을 후처리로 붙여, 충격파와 함께 있는 계면을 상수 두께로 유지한다.

왜 계면은 저절로 뚱뚱해지는가#

5-방정식 모델은 보존형에 색함수(color function) 하나를 얹은 준보존 시스템이다. 부피분율 $\alpha_1$ (셀 안에서 유체 1이 차지하는 비율)의 이송 방정식은 이렇게 생겼다.

\frac{\partial \alpha_1}{\partial t} + \mathbf{u}\cdot\nabla \alpha_1 = \alpha_1 K\, \nabla\cdot\mathbf{u}

$\alpha_1$ 은 유체 1의 부피분율, $\mathbf{u}$ 는 혼합 속도, 우변의 Kapila 항은 두 상의 압축률 차이를 보정한다.

문제는 이 방정식을 충격파 포착 스킴으로 풀 때 생긴다. 셀 경계에서 값을 매끄럽게 보간하는 순간, 날카로워야 할 $\alpha_1$ 의 계단이 매 스텝 조금씩 번진다. 이 수치 확산은 누적된다. 스텝이 쌓이면 계면은 몇 칸에서 수십 칸으로 퍼지고, 기포 모양·충격파 반사 위치 같은 핵심 특징이 통째로 사라진다.

WENO 같은 고차 스킴은 확산을 줄이지만 진동을 낳고, 다차원에서는 비싸다. 그래서 저자들은 방향을 튼다. 스킴 자체를 건드리지 않고, 매 스텝이 끝난 뒤 $\alpha$ 장을 따로 손질하기로 한다.

안티확산: 샤프닝 방정식#

쓰는 식은 Shukla 등(2010)의 정규화 방정식이다. 실제 물리 시간이 아닌 가상 시간 $\tau$ 에 대해 $\alpha$ 를 몇 번 진화시킨다.

\frac{\partial \alpha}{\partial \tau} = \nabla\cdot\!\left(\varepsilon\, \nabla \alpha\right) - \nabla\cdot\!\left(a\, \alpha(1-\alpha)\, \hat{\mathbf{n}}\right)

$\varepsilon$ 은 두께를 규제하는 확산 계수, $a$ 는 계면을 조이는 압축 세기, $\hat{\mathbf{n}} = \nabla\alpha / |\nabla\alpha|$ 는 계면 법선 방향이다.

두 항이 줄다리기를 한다. 첫 항은 보통의 확산이라 $\alpha$ 를 퍼뜨린다. 둘째 항은 부호가 반대인 안티확산(anti-diffusion) 이라 계면을 법선 방향으로 조인다. $\alpha(1-\alpha)$ 인자 덕에 압축은 계면( $0<\alpha<1$ )에서만 작동하고, 순수 유체 영역( $\alpha=0$ 또는 $1$ )에서는 꺼진다.

정상상태( $\partial_\tau \alpha = 0$ )를 1차원에서 풀면 두 플럭스가 균형을 이룬다.

\varepsilon\, \frac{d\alpha}{dx} = a\, \alpha(1-\alpha) \quad\Longrightarrow\quad \alpha(x) = \tfrac{1}{2}\!\left[1 + \tanh\!\left(\frac{a}{2\varepsilon}\,x\right)\right]

계면은 하이퍼볼릭 탄젠트 모양으로 수렴하고, 그 두께는 $\varepsilon/a$ 에 비례한다. 즉 두께를 사용자가 정한다. 이게 저자들이 강조하는 "상수 두께" 성질이다. 계면이 항상 같은 칸 수를 유지하면 충격·접촉 불연속을 안정적으로 포착할 수 있다.

1D로 직접: 번짐과 샤프닝의 줄다리기#

위 화면의 줄다리기를 같은 식으로 적으면 아래다. 이송이 남긴 번짐을 흉내 낸 뒤, 샤프닝을 가상 시간으로 몇 번 돌린다.

import numpy as np
 
N, dx = 201, 1.0 / 200
x = np.linspace(0, 1, N)
 
def band_width(a, lo=0.05, hi=0.95):
    # 계면으로 간주하는 셀 수 (두께 측정)
    return int(np.sum((a > lo) & (a < hi)))
 
def smear_once(a, D=0.16):
    # 이송 스킴이 매 스텝 계면을 번지게 하는 효과
    lap = np.zeros_like(a)
    lap[1:-1] = a[2:] - 2 * a[1:-1] + a[:-2]
    return a + D * lap
 
def sharpen_sweep(a, comp_a, eps, iters=200):
    # 플럭스형 안티확산 정규화 (Shukla 2010, 논문 식 (38)-(39))
    dtau = 0.9 * min(dx * dx / (2 * eps), dx / comp_a)
    for _ in range(iters):
        af = 0.5 * (a[:-1] + a[1:])          # 면 중심 alpha
        dA = a[1:] - a[:-1]
        s = np.sign(dA)
        J = comp_a * af * (1 - af) * s - eps * dA / dx   # 면 플럭스
        a[1:-1] = np.clip(a[1:-1] - dtau / dx * (J[1:] - J[:-1]), 0, 1)
    return a
 
# 날카로운 계단에서 시작 -> 20번 번지게 둔다
a = 0.5 * (1 + np.tanh((x - 0.5) / 0.012))
for _ in range(20):
    a = smear_once(a)
print("번진 뒤 밴드:", band_width(a), "칸")     # -> 넓어짐
 
a = sharpen_sweep(a, comp_a=1.0, eps=0.004)
print("샤프닝 뒤 밴드:", band_width(a), "칸")   # -> 목표 두께로 수축

출력은 번진 뒤 밴드가 넓고, 샤프닝 뒤에는 $\varepsilon/a$ 가 정한 몇 칸으로 줄어든다. 스킴은 그대로다. sharpen_sweep은 결과 배열만 받는다. 리만 솔버가 무엇이었는지는 모른다.

기본 인자는 comp_a=1.0, eps=0.004다. 두께 척도 $\varepsilon/a$ 는 0.004. 정상 프로파일의 기울기 계수는 $a/(2\varepsilon)=125$ 다. 코드가 시작한 tanh의 분모는 0.012라, 목표가 초기 계단보다 더 가파르다. dx는 $1/200=0.005$ . 가상시간 상한은 $dx^2/(2\varepsilon)=0.000025/0.008=0.003125$ 와 $dx/a=0.005$ 중 작은 쪽이다. 확산 제한이 이긴다. 0.9를 곱한 $d\tau$ 는 0.0028125. 이 상한을 넘기면 안티확산이 부호를 뒤집으며 진동한다.

면 플럭스의 압축 조각은 $\alpha(1-\alpha)$ 에 비례한다. 그 곱은 $\alpha=0.5$ 에서 0.25로 최대고, 순수 유체에서는 0이라 압축이 꺼진다. 밴드 함수는 0.05와 0.95 사이를 센다. $\varepsilon$ 만 두 배면 목표 두께도 두 배, $a$ 만 두 배면 절반이다. 도입의 200스텝은 물리 시간이다. 6칸이 20칸이 되면 밴드는 약 3.3배다. smear_once의 D=0.16과 샤프닝의 iters=200은 다른 층이다. 논문의 스텝당 1~3회는 물리 스텝의 확산이 작을 때의 이야기다. Kapila 우변 $\alpha_1 K\nabla\cdot\mathbf{u}$ 는 충격에서 발산이 스파이크일 때 부피분율을 직접 민다. 수치 확산과는 별개의 소스라, IST는 그 소스를 지우지 않는다. 5-방정식의 압력 이완을 따로 두는 Saurel(2009) 절차와 달리, 여기서는 계면 칸 수만 되돌린다.

두께를 상수로 만든다는 것#

단순히 계면을 얇게 만드는 것과, 두께를 상수로 고정하는 것은 다르다. 저자들이 이전 접근(Tiwari 등 2013)과 선을 긋는 지점이 여기다. 소스 항에 샤프닝 함수를 더하는 방식은 확산 오차를 줄이지만, 시간·공간에 따라 두께가 들쭉날쭉해진다. 반면 후처리로 $\varepsilon/a$ 를 고정하면 어디서나 같은 두께가 나온다.

압축성 유동에서는 여기에 한 겹이 더 붙는다. $\alpha$ 를 손질하면 혼합 밀도·에너지가 함께 바뀌어야 열역학이 일관된다. 저자들은 $\alpha$ 만 고치는 게 아니라, 보존 변수 전체( $\alpha_1\rho_1$ , $\alpha_2\rho_2$ , 운동량, 에너지)를 혼합 규칙에 맞게 재분배한다. 이 "혼합 일관 정규화" 덕에, 계면에서 속도·압력·온도 평형이 깨지지 않는다. 논문의 1D 순수 계면 이송 테스트가 이걸 보여준다. 물과 공기가 같은 속도로 흐를 때, 샤프닝을 켜도 압력·온도가 그대로 유지된다.

2D 기포로: 원형 계면 유지#

1D에서 통했으니 2D 기포로 넘어간다. 논문의 shock-bubble 그림은 계면이 몇 칸 두께의 원을 끝까지 유지한다. numerical diffusion을 켜면 테두리가 흐려지고, sharpening과 compression a를 올리면 다시 링으로 조여진다.

compression a: 1.20

interface band: 0 cells

teal ring = cells with 0.1 < α < 0.9. Fewer = sharper.

청록 링이 조여도 기포 중심은 그대로다. 커지거나 줄지 않는다. 줄어드는 것은 링의 칸 수뿐이다. 압축이 $\hat{\mathbf{n}}$ 을 따라가서 곡률이 찌그러지지 않는다.

재현하며 부딪힌 것#

$|\nabla\alpha|\to 0$ 인 순수 유체에서 법선은 0으로 나눈다. 논문은 $a\,\alpha(1-\alpha)$ 가 거기서 0이라 넘어가지만, 분모에 작은 수를 더하지 않으면 NaN이 난다. $d\tau$ 가 크면 계면이 진동한다. CFL형 상한이 필요하다. 논문의 스텝당 1~3회는 물리 스텝의 확산이 작을 때의 이야기다. 성긴 격자나 강한 충격은 그 횟수로 부족할 수 있다.

한 가지 더. 후처리라는 우아함에는 대가가 있다. 샤프닝은 지배 방정식과 별개로 $\alpha$ 를 바꾸므로, 그 순간의 이송은 원래 PDE의 해가 아니다. 두께가 상수로 유지되는 대신, 계면 근처의 국소 보존이 미세하게 어긋날 여지가 있다. 저자들은 혼합 일관 재분배로 이를 최소화하지만 완전히 0은 아니다. OpenFOAM의 interfoam이 쓰는 압축 항( $\nabla\cdot(\alpha(1-\alpha)\mathbf{u}_r)$ )도 같은 계열의 아이디어라, 이 기법을 압축성 솔버에 이식할 때 비교 기준으로 삼을 만하다.

ε/a를 잠그지 않으면 칸 수가 다시 뜬다#

소스항에 샤프닝을 넣으면 확산 오차는 줄어도 두께가 시공간에서 흔들린다. 스텝이 끝난 뒤 $\varepsilon/a$ 를 고정해야 충격 옆에서도 같은 칸 수가 남는다. $\alpha$ 만 고치고 보존 변수를 혼합 규칙으로 다시 나누지 않으면, 그 칸에서 압력과 온도 평형이 먼저 나간다.

$d\tau$ 가 $dx^2/(2\varepsilon)$ 보다 크면 안티확산은 계면을 조이는 대신 진동시킨다. 두께를 상수로 사는 대가는, 그 순간 이송이 원래 PDE의 해가 아니라는 점이다. 국소 보존이 미세하게 어긋나는 자리는 Abgrall 기준이 잡는 비보존항과 같은 종류다.