격자를 절반으로 줄였더니 발산이 두 배 빨리 왔다 — 6방정식 two-fluid 모델의 복소 고유값

복소 고유값은 격자를 조밀하게 할수록 더 빨리 터진다. 고칠 자리는 이산화가 아니라 계면압력 닫힘 관계다.

격자를 절반으로 줄였더니 발산이 두 배 빨리 왔다#

수치 코드에서 답이 터지면 보통 격자를 조밀하게 해본다. 오차가 줄면 이산화 문제, 그대로면 물리 모델 문제다. 대략 이 순서로 좁혀 들어간다.

그런데 반대로 움직이는 경우가 있다. 격자를 절반으로 줄였더니 발산이 정확히 두 배 빨리 왔다. 셀 수를 네 배로 늘리면 네 배 빨리 온다. 시간 스텝을 줄여도 성장률은 그대로다.

이 증상은 이산화 버그가 아니다. 지배방정식 자체가 초기값 문제로서 부적절(ill-posed, 초기 교란이 파장에 반비례해 무한히 빠르게 자라는 상태)하다는 신호다. Pandare와 Luo가 2018년 AIAA 논문에서 밀도기반 유한체적 two-fluid 솔버를 세울 때 가장 먼저 손댄 것도 이 자리였다.

이 글은 단일압력 6방정식 two-fluid 모델의 고유값을 직접 뽑아, 어디서 복소수가 생기고 계면압력 항의 계수가 얼마여야 실수로 돌아오는지를 숫자로 확인한다. 답은 정확히 1이다.

느린 파동 한 쌍이 실수축을 떠난다#

두 상을 서로 침투하는 연속체로 보는 것이 two-fluid 모델이다. 상별로 질량·운동량·에너지를 따로 푼다. 두 상의 압력을 하나로 묶으면( $p_g = p_l \equiv p$ ) 6개 PDE가 남는다. Wallis 모델, 또는 단일압력 6방정식 모델이라 부른다.

1차원에서 압축성을 잠시 끄고 원시변수 $(\alpha_g, p, u_g, u_l)$ 로 준선형화하면, 느린 파동 한 쌍의 고유값이 닫힌 형태로 나온다.

\lambda_\pm = \frac{\alpha_l \rho_g u_g + \alpha_g \rho_l u_l}{\alpha_l \rho_g + \alpha_g \rho_l} \pm \frac{\sqrt{(\sigma - 1)\, \alpha_g \alpha_l \rho_g \rho_l}\; |u_l - u_g|}{\alpha_l \rho_g + \alpha_g \rho_l}

$\alpha_k$ 는 체적분율, $\rho_k$ 는 상밀도, $u_k$ 는 상속도다. $\sigma$ 는 곧 설명할 계면압력 항의 계수다. 앞의 항은 밀도로 가중한 평균 속도, 뒤의 항이 두 파동이 갈라지는 폭이다.

근호 안이 전부다. $\sigma < 1$ 이면 음수가 되고, 두 고유값은 켤레복소수 한 쌍이 된다. 미끄럼 속도 $|u_l - u_g|$ 가 0이 아닌 한 이 일은 반드시 일어난다. 즉 두 상이 서로 다른 속도로 흐르는 순간 모델이 부적절해진다.

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

sigma0.00slip u_r10.0 m/salpha_g0.50

Push sigma up from 0. The two red dots slide down the imaginary axis, meet at σ = 1, then split along the real axis and turn green — and the right panel stops growing (envelope now ×1.00) and starts propagating as two separate void waves. Set slip u_r to 0 and the whole pair collapses onto one point: no slip, no problem.

sigma를 0에서 올리면 왼쪽 복소평면의 붉은 점 두 개가 허수축을 따라 내려와 1에서 만나고, 그다음 실수축 위로 갈라지며 초록으로 바뀐다. 오른쪽 교란이 자라기를 멈추고 좌우로 흐르기 시작하는 순간이 같은 순간이다. slip u_r을 0으로 내리면 문제 자체가 사라지는 것도 확인해보자.

표 한 장 — 7방정식·6방정식·계면압력의 세 칸#

이 자리를 두고 세 가지 선택지가 있다. 세로로 세워 비교하면 각각이 무엇을 사고 무엇을 파는지가 드러난다.

	7방정식 (Baer–Nunziato)	6방정식 원형 (Wallis)	6방정식 + 계면압력
압력	상별로 둘	하나	하나
고유값	항상 실수	미끄럼 있으면 복소	$\sigma \ge 1$ 이면 실수
미지수	상분율 이송식 추가	최소	최소
대가	압력 완화 항, 강성	부적절	물리적 근거가 약한 $\sigma$
적용 한계	조밀 입자·현탁액에서 물리적	그대로는 못 씀	공학적 절충

7방정식 모델은 체적분율에 자체 이송 방정식을 하나 더 준다. 그 대가로 쌍곡성을 확보하지만, 압력 완화 항이 강성을 들여온다. flux splitting으로 Baer–Nunziato를 다뤘던 글에서 그 구조를 정리한 적이 있다. 문제는 이 모델이 물리적으로 정당한 범위가 조밀하게 채워진 입자·현탁액 쪽이라는 점이다. 물과 공기가 층을 이루고 흐르는 배관에는 잘 맞지 않는다.

Python으로 뽑은 4×4 고유값#

닫힌 형태를 믿기 전에 원래 시스템을 그대로 풀어본다. 압축성을 살려 $A W_t + B W_x = 0$ 을 세우고 $A^{-1}B$ 의 고유값을 뽑는다. 물과 공기, $\alpha_g = 0.5$ , 기체가 10 m/s로 앞선다.

import numpy as np
 
def interfacial_dp(a, rg, rl, ur, sigma):
    """Stuhmiller 보정: p_int = p - dp"""
    return sigma * a * (1 - a) * rg * rl * ur**2 / (a * rl + (1 - a) * rg)
 
def two_fluid_matrices(a, rg, rl, cg, cl, ug, ul, sigma):
    """A W_t + B W_x = 0,  W = (alpha_g, p, u_g, u_l)"""
    dp = interfacial_dp(a, rg, rl, ul - ug, sigma)
    kg, kl = a / (rg * cg**2), (1 - a) / (rl * cl**2)
    A = np.array([[ 1.0, kg,  0.0,    0.0],
                  [-1.0, kl,  0.0,    0.0],
                  [ 0.0, 0.0, a * rg, 0.0],
                  [ 0.0, 0.0, 0.0,    (1 - a) * rl]])
    B = np.array([[ ug,  ug * kg, a,          0.0],
                  [-ul,  ul * kl, 0.0,        1 - a],
                  [ dp,  a,       a * rg * ug, 0.0],
                  [-dp,  1 - a,   0.0,        (1 - a) * rl * ul]])
    return A, B
 
def char_speeds(sigma, a=0.5, rg=1.2, rl=1000.0, cg=340.0, cl=1500.0, ug=10.0, ul=0.0):
    A, B = two_fluid_matrices(a, rg, rl, cg, cl, ug, ul, sigma)
    return np.linalg.eigvals(np.linalg.solve(A, B))
 
print("air/water, alpha_g=0.5, u_g=10, u_l=0 m/s")
print("sigma   max|Im lambda|   slow pair Re")
for s in [0.0, 0.5, 0.9, 1.0, 1.1, 1.5]:
    lam = char_speeds(s)
    slow = np.sort(lam.real)[1:3]
    print("%5.2f   %12.5f   %8.4f %8.4f" % (s, np.abs(lam.imag).max(), slow[0], slow[1]))

air/water, alpha_g=0.5, u_g=10, u_l=0 m/s
sigma   max|Im lambda|   slow pair Re
 0.00        0.34614     0.0120   0.0120
 0.50        0.24481     0.0120   0.0120
 0.90        0.10967     0.0120   0.0120
 1.00        0.00718     0.0120   0.0120
 1.10        0.00000    -0.0972   0.1212
 1.50        0.00000    -0.2326   0.2566

$\sigma = 0$ 에서 허수부가 0.346 m/s다. 두 느린 파동의 실수부는 0.0120으로 붙어 있다. 기체가 10 m/s로 흐르는데 파동 속도가 0.012 m/s인 이유는 밀도 가중 때문이다. 물이 공기보다 830배 무거우니 평균이 물 쪽으로 끌려간다.

$\sigma$ 를 올리면 허수부가 줄고, 1.1에서 0이 되면서 두 파동이 $-0.097$ 과 $0.121$ 로 갈라진다. $\sigma = 1.0$ 에서 아직 0.00718이 남은 것은 압축성 때문이다. 닫힌 형태는 비압축 극한에서 유도했으므로, 유한한 음속이 임계값을 1보다 아주 조금 위로 밀어 올린다.

닫힌 형태가 말하는 임계값은 정확히 1이다#

이제 닫힌 형태로 같은 값을 재보고, 임계 $\sigma$ 를 이분법으로 찾는다.

from math import sqrt, pi
 
def material_pair(a, sigma, rg=1.2, rl=1000.0, ug=10.0, ul=0.0):
    """느린(물질) 파동 한 쌍의 비압축 극한 닫힌 형태"""
    al = 1.0 - a
    den = al * rg + a * rl
    mean = (al * rg * ug + a * rl * ul) / den
    disc = (sigma - 1.0) * a * al * rg * rl * (ul - ug) ** 2 / den**2
    if disc >= 0.0:
        return (mean - sqrt(disc), mean + sqrt(disc)), 0.0
    return (mean, mean), sqrt(-disc)
 
print("closed form vs the 4x4 eigenvalues above")
for s in [0.0, 0.5, 0.9, 1.1, 1.5]:
    (r1, r2), im = material_pair(0.5, s)
    print("sigma=%4.2f  Re = %8.4f %8.4f   |Im| = %8.5f" % (s, r1, r2, im))
 
print()
print("growth rate of the shortest resolved mode, L = 1 m, sigma = 0")
_, im0 = material_pair(0.5, 0.0)
for n in [50, 100, 200, 400, 800]:
    k = pi * n          # k = pi / dx, dx = 1/n
    print("N=%4d  dx=%7.5f  k=%8.1f 1/m  growth=%8.2f 1/s" % (n, 1.0 / n, k, k * im0))
 
print()
print("critical sigma (incompressible limit) for a few states")
for a in [0.1, 0.5, 0.9]:
    for ur in [1.0, 30.0]:
        lo, hi = 0.0, 5.0
        for _ in range(60):
            mid = 0.5 * (lo + hi)
            _, im = material_pair(a, mid, ug=ur)
            if im > 0.0: lo = mid
            else: hi = mid
        print("alpha_g=%.1f  u_r=%4.1f  ->  sigma_c = %.6f" % (a, ur, hi))

closed form vs the 4x4 eigenvalues above
sigma=0.00  Re =   0.0120   0.0120   |Im| =  0.34599
sigma=0.50  Re =   0.0120   0.0120   |Im| =  0.24466
sigma=0.90  Re =   0.0120   0.0120   |Im| =  0.10941
sigma=1.10  Re =  -0.0974   0.1214   |Im| =  0.00000
sigma=1.50  Re =  -0.2327   0.2566   |Im| =  0.00000
 
growth rate of the shortest resolved mode, L = 1 m, sigma = 0
N=  50  dx=0.02000  k=   157.1 1/m  growth=   54.35 1/s
N= 100  dx=0.01000  k=   314.2 1/m  growth=  108.70 1/s
N= 200  dx=0.00500  k=   628.3 1/m  growth=  217.40 1/s
N= 400  dx=0.00250  k=  1256.6 1/m  growth=  434.79 1/s
N= 800  dx=0.00125  k=  2513.3 1/m  growth=  869.58 1/s
 
critical sigma (incompressible limit) for a few states
alpha_g=0.1  u_r= 1.0  ->  sigma_c = 1.000000
alpha_g=0.1  u_r=30.0  ->  sigma_c = 1.000000
alpha_g=0.5  u_r= 1.0  ->  sigma_c = 1.000000
alpha_g=0.5  u_r=30.0  ->  sigma_c = 1.000000
alpha_g=0.9  u_r= 1.0  ->  sigma_c = 1.000000
alpha_g=0.9  u_r=30.0  ->  sigma_c = 1.000000

닫힌 형태와 4×4 고유값이 소수 셋째 자리까지 일치한다. 체적분율을 0.1에서 0.9로 바꾸고 미끄럼을 1에서 30 m/s로 바꿔도 임계값은 1.000000이다. Stuhmiller가 제안한 보정

p_{\text{int}} = p - \sigma\, \frac{\alpha_g \alpha_l \rho_g \rho_l}{\alpha_g \rho_l + \alpha_l \rho_g}\, u_r^2

에서 $\sigma = 1$ 이 임의의 튜닝 값이 아닌 이유가 여기 있다. 근호 안을 정확히 0으로 만드는 최소 계수다. 실무에서는 여유를 두고 1보다 조금 큰 값을 쓴다.

불안정과 부적절 문제는 다르다#

수치적으로 불안정한 스킴은 시간 스텝을 줄이면 나아진다. 부적절 문제는 그렇지 않다. 성장률이 파수에 비례하기 때문이다.

\text{growth}(k) = k \, |\mathrm{Im}\,\lambda|, \qquad k = \frac{\pi}{\Delta x}

$\Delta x$ 를 절반으로 줄이면 표현 가능한 최단 파장이 절반이 되고, 성장률이 두 배가 된다. 위 출력에서 $N = 50$ 의 54.35 1/s가 $N = 800$ 에서 869.58 1/s로 정확히 16배가 된 것이 그것이다. 격자를 조밀하게 할수록 답이 더 빨리 죽는다.

sigma0.00slip u_r10.0 m/st = 0.0 ms

Watch the order in which the lanes hit the blow-up line: the finest grid always gets there first, and doubling N halves the time. Drag sigma past 1 and every lane goes flat at the same instant — the cure is in the closure, not in the mesh.

네 개의 격자가 같은 교란을 안고 동시에 출발한다. 어느 레인이 먼저 blow-up 선에 닿는지, 그리고 sigma를 1 위로 올렸을 때 네 레인이 동시에 평평해지는지를 보자. 고치는 자리가 격자가 아니라 닫힘 관계임이 그림 하나로 드러난다.

실제 코드에서 이 증상이 가려지는 경우도 많다. 1차 상류 차분의 수치 확산이 $O(k^2 \Delta x)$ 로 감쇠를 넣어주면, 성장률이 상쇄되어 계산이 그럭저럭 굴러간다. 그래서 저차 스킴에서 멀쩡하던 코드가 고차로 올리는 순간 터진다. 보존형과 원시형이 갈라지던 자리와 비슷한 구조다. 수치 확산이 모델의 결함을 대신 갚고 있었을 뿐이다.

표의 나머지 칸 — 밀도기반이 낮은 마하수에서 살아남는 법#

쌍곡성을 되찾았다고 끝이 아니다. 다상 유동의 실제 응용은 대부분 마하수가 매우 낮다. 밀도기반 솔버는 이 영역에서 음속 CFL에 묶여 시간 스텝이 무너진다.

전통적으로 이 자리는 압력기반 방법이 차지해왔다. 속도장을 솔레노이드로 가정해 음속을 방정식에서 지워버리므로 CFL이 유속만으로 정해진다. 대신 압축성을 엄밀하게 다루지 못한다. 비등처럼 고온 현상이 들어오면 오차가 커진다.

Pandare와 Luo가 택한 길은 밀도기반을 유지하되 원시변수 $[p, v, T]$ 로 변환해 완전 음해로 푸는 것이다. 압력을 미지수로 세우면 낮은 마하수에서 조건수가 좋아진다. 항력·가상질량 같은 계면 힘 항도 음해로 처리해 시간 스텝 제약을 더 푼다.

플럭스 쪽에도 같은 절충이 있다. 강한 충격파가 물질 계면과 만나면 AUSM $^+$ -up이 음압을 낸다. 기존 해법은 그 자리에만 정확 리만 해법기를 부르는 것이었지만, 뉴턴 반복 비용이 크다. 논문은 질량 플럭스에 체적분율 결합 항을 하나 넣어 같은 강건성을 얻는다. Lax–Friedrichs 성 소산을 체적분율 점프에 비례해 넣는 셈이다. 정지한 계면을 흔들지 않아야 한다는 조건은 계면 포착 스킴에서 CFL 상한을 재던 글에서도 같은 이름으로 등장했다.

세 칸 중 어디에 서 있는지부터 확인한다#

two-fluid 솔버를 새로 켤 때, 격자나 스킴을 만지기 전에 확인할 것이 셋 있다.

첫째, 미끄럼 속도가 0이 아닌 상태에서 자코비안 고유값을 뽑아본다. 4×4 행렬 하나면 된다. 허수부가 나오면 이산화로 해결할 문제가 아니다.

둘째, 격자를 두 배로 조밀하게 하고 발산 시점을 잰다. 시점이 절반이 되면 부적절, 늦어지면 이산화 문제다. 이 한 번의 실험이 진단을 가른다.

셋째, 계면압력 계수를 코드에서 찾아 값을 확인한다. 1보다 작으면 그 코드는 수치 확산으로 버티고 있는 것이다. 고차로 올리기 전에 이 값을 먼저 올려야 한다.

격자를 절반으로 줄였더니 발산이 두 배 빨리 왔다 — 6방정식 two-fluid 모델의 복소 고유값