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

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

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

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

sigma를 0에서 1.1로 밀면 왼쪽 복소평면의 점이 실수축으로 내려앉고, slip u_r을 0으로 내리면 교란이 자라지 않는다.

점이 초록으로 바뀌는 순간이 오른쪽 교란이 멈추는 순간이다. 조밀한 격자의 레인이 더 먼저 천장에 닿는지도 같은 화면에서 보인다.

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

단일압력 6방정식의 고유값을 직접 뽑아 봤다. 복소수가 생기는 자리와, 계면압력 계수가 얼마일 때 실수축으로 돌아오는지를 숫자로 맞춘다. 답은 정확히 1이다.

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

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

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

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

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

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

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

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보다 아주 조금 위로 밀어 올린다.

밀도 가중의 나눗셈을 다시 적으면 0.0120이 어디서 왔는지 보인다. $\alpha_g = 0.5$, $\rho_g = 1.2$, $\rho_l = 1000$, $u_g = 10$이면 분모는 $0.5\times 1.2 + 0.5\times 1000 = 500.6$이다. 분자는 $0.5\times 1.2\times 10 = 6$이라 평균 속도는 $6/500.6 = 0.0120$이다. 물이 공기보다 $1000/1.2 \approx 833$배 무겁다는 말이 이 한 줄이다. 기체가 10 m/s로 달려도 느린 파동의 실수부는 0.012 m/s에 붙어 있다.

허수부 비율도 닫힌 형태 그대로다. $\sigma = 0$의 0.34614와 $\sigma = 0.5$의 0.24481의 비는 $\sqrt{2} \approx 1.414$다. $\sigma = 0.9$에서는 $\sqrt{0.1}$만큼 줄어 0.10967이 된다. $|\mathrm{Im}\,\lambda|$가 $\sqrt{1-\sigma}$에 비례한다. $\sigma = 1.1$에서 실수부가 $-0.0972$와 $0.1212$로 갈라져도 중점은 0.012다. 반폭 0.109는 $\sigma = 0.9$의 허수부와 같은 크기다. 근호 안의 $(\sigma-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

sigma를 0에 두면 셀이 많은 레인이 먼저 blow-up 선에 닿고, 1을 넘기면 네 레인이 같이 눕는다. N = 50의 54.35 1/s와 N = 800의 869.58 1/s가 같은 출발에서 갈린다.

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

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

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

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

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

표의 곱셈도 같은 비례다. $N = 50$에서 $k = 157.1$ 1/m이고 $|\mathrm{Im}\,\lambda| = 0.34599$이면 $157.1 \times 0.34599 \approx 54.35$ 1/s다. $N$을 두 배 할 때마다 $k$도 두 배라 54.35, 108.70, 217.40, 434.79, 869.58로 정확히 두 배씩 간다. $N = 50$에서 $N = 800$은 격자 16배, 성장률 16배다. $\Delta x = 0.02$ m일 때 초당 54번 불어나던 교란은 $\Delta x = 0.00125$ m이면 초당 870번에 가깝다. 시간 스텝을 줄여도 파수 157.1은 그대로다. 줄어드는 것은 $\Delta t$이지 $k$가 아니다.

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

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

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

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

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

플럭스 쪽에도 같은 절충이 있다. 강한 충격파가 물질 계면과 만나면 AUSM$^+$-up이 음압을 낸다. 기존 해법은 그 자리에만 정확 리만 해법기를 부르는 것이었지만, 뉴턴 반복 비용이 크다. 논문은 질량 플럭스에 체적분율 결합 항을 하나 넣어 같은 강건성을 얻는다. Lax–Friedrichs 성 소산을 체적분율 점프에 비례해 넣는 셈이다. 정지한 계면을 흔들지 않아야 한다는 조건은 번진 계면을 매 스텝 되돌리는 안티확산과도 같은 제약이다.

σ가 1 아래면 차수를 올리지 않는다#

미끄럼이 남은 상태에서 4×4를 한 번 풀어 허수부가 보이면, 격자를 반으로 줄이는 실험은 진단용일 뿐이다. $N = 50$의 54.35 1/s가 $N = 800$에서 869.58 1/s가 되는 계산은 스킴이 아니라 모델이 병이다. 코드에서 $\sigma$를 먼저 읽는다.

1보다 작으면 그 런은 1차 상류의 $O(k^2 \Delta x)$ 수치 확산으로 버티는 중이다. 차수를 올리면 그 완충이 빠진다. $\sigma$를 1.1 근처로 올린 뒤에야 격자 조밀화가 오차를 줄이는 쪽으로 돌아선다. 저마하에서 밀도기반으로 버틴 다음의 자리는 Poisson 없이 압력을 세우는 명시적 결합을 따로 저울질할 일이다.