K를 지웠더니 물이 1.86배 압축됐다 — 격자를 조여도 줄지 않는 계면 오차

계면 셀에서 ψ+K는 거의 0이다. K는 ψ∇·u를 지우려고 들어간 계수이고, 그 항을 빼면 액체가 기체의 압축률로 눌린다.

압축성 VOF의 색상 함수 이송식에는 항이 하나 더 붙는다. $-(\psi+K)\nabla\cdot u$ 다. 그 $K$ 를 코드에 옮기고 값을 찍어봤다. 물과 공기의 계면 셀에서 $\psi=0.5$ 일 때 $K=-0.499947$ , 그래서 $\psi+K=5.3\times10^{-5}$ 였다. 식을 정확하게 만들려고 넣은 계수가 그 항을 거의 0으로 지워버린 것이다.

이 글은 Denner의 압력기반 완전연성 압축성 다상 유한체적법 논문에서 $K$ 를 등압 종결로 직접 유도한다. 그리고 $K$ 를 0으로 두거나 부호를 뒤집었을 때 1차원 압축 문제에서 액체 체적이 각각 1.86배·4.23배 더 줄어드는 것을 재현한다. 마지막으로 그 오차가 격자를 8배 조여도 7%밖에 줄지 않는 이유를 본다.

압축성 VOF에는 왜 $-(\psi+K)\nabla\cdot u$ 가 붙어 있나#

비압축성 VOF에서 색상 함수는 그냥 흘러간다. $\nabla\cdot u=0$ 이므로 이송만 하면 끝이다. 압축성이 되면 유체가 늘어나고 줄어든다. 색상 함수는 체적률이므로 그 팽창·수축을 따라가야 한다.

논문의 식은 이렇다.

\frac{\partial \psi}{\partial t} + \frac{\partial u_i \psi}{\partial x_i} - (\psi + K)\frac{\partial u_i}{\partial x_i} = 0

$\psi$ 는 액체의 체적률, $u_i$ 는 속도, $K$ 는 물질에 따라 달라지는 압축률 인자(compressibility factor)다. 세 번째 항을 옮기면 물질미분 형태가 나온다.

\frac{\mathrm{D}\psi}{\mathrm{D}t} = K\,\nabla\cdot u

$K=0$ 이면 체적률은 물질을 따라 변하지 않는다. 그럴듯해 보이지만 틀렸다. 계면 셀의 $\nabla\cdot u$ 는 혼합물 전체의 값이고, 혼합물의 압축률은 부드러운 쪽이 지배한다.

숫자로 보자. 1기압 물의 체적탄성계수는 $\rho_l a_l^2 = 2.64\times10^9$ Pa, 공기는 $\rho_g a_g^2 = 1.40\times10^5$ Pa다. 비는 18,857배다. 혼합물의 압축률은 체적률 가중 합이므로, $\psi=0.5$ 에서도 공기 항이 거의 전부를 차지한다. 그 결과 혼합물의 음속은 23.65 m/s까지 떨어진다. 물보다 69배, 공기보다 14배 느리다.

아래 시뮬레이션에서 체적률을 직접 움직여보자.

ψ (liquid volume fraction) = 0.500pressure = 1 bar

a_mix = 23.65 m/s  |  ρ_mix = 500.60 kg/m³  |  a_mix / a_gas = 0.069
crossing 2.0 m: liquid 1.23 ms  |  gas 5.85 ms  |  mixture 84.57 ms

$\psi$ 슬라이더를 0 근처에서 조금만 올려보면 음속이 341 m/s에서 100 m/s 아래로 수직 낙하한다. 아래 트랙은 같은 2 m 구간을 음향 펄스가 건너는 데 걸리는 시간이다. 액체 비율 0.1%만 섞여도 혼합물 레인이 기체 레인보다 눈에 띄게 뒤처진다.

$K$ 는 어디서 나오는 숫자인가#

$K$ 는 가정에서 나온다. 등압 종결(isobaric closure) — 계면 셀 안의 두 상이 같은 압력을 공유하고, 각자 자기 체적탄성계수로 반응한다는 가정이다. 논문이 ACID에서 쓰는 국소 단상 가정과 같은 내용이다.

라그랑주 셀을 따라가며 상 $k$ 의 질량이 보존된다고 쓰면

\frac{1}{V_k}\frac{\mathrm{D}V_k}{\mathrm{D}t} = -\frac{1}{\rho_k a_k^2}\frac{\mathrm{D}p}{\mathrm{D}t}

가 된다. $V_k$ 는 그 상이 차지하는 체적, $\rho_k a_k^2$ 는 체적탄성계수다. 두 상을 더하면 혼합물의 발산이 나온다.

\nabla\cdot u = -\left[\frac{\psi}{\rho_l a_l^2} + \frac{1-\psi}{\rho_g a_g^2}\right]\frac{\mathrm{D}p}{\mathrm{D}t}

대괄호 안이 혼합물의 압축률이다. 여기서 Wood 음속이 자동으로 나온다. 이제 $\psi = V_l/V$ 를 미분하고 위 두 식을 넣으면

\frac{\mathrm{D}\psi}{\mathrm{D}t} = \psi(1-\psi)\left[\frac{1}{\rho_g a_g^2} - \frac{1}{\rho_l a_l^2}\right]\frac{\mathrm{D}p}{\mathrm{D}t}

가 남는다. 두 압축률의 차이가 곱해져 있다. 두 상의 체적탄성계수가 같으면 체적률은 변하지 않는다. $\mathrm{D}\psi/\mathrm{D}t = K\nabla\cdot u$ 와 비교해 $\mathrm{D}p/\mathrm{D}t$ 를 소거하면 $K$ 가 떨어진다.

K = \frac{\rho_g a_g^2 - \rho_l a_l^2}{\dfrac{\rho_g a_g^2}{1-\psi} + \dfrac{\rho_l a_l^2}{\psi}}

분모에 $1/\psi$ 와 $1/(1-\psi)$ 가 있다. 순수 셀에서 발산하므로 코드에는 약분한 형태로 적어야 한다.

K = \frac{(\rho_g a_g^2 - \rho_l a_l^2)\,\psi(1-\psi)}{\rho_g a_g^2\,\psi + \rho_l a_l^2\,(1-\psi)}

이 형태는 $\psi=0$ 과 $\psi=1$ 에서 그냥 0이 된다. 분자의 부호는 $\psi$ 가 어느 상을 가리키는지에 전부 달려 있다. 논문 식을 그대로 베끼면 라벨이 바뀌어 부호가 뒤집힐 수 있다. 뒤에서 그 결과를 본다.

$\psi+K$ 가 $5\times10^{-5}$ 라면 두 항을 지워도 되나#

$\rho_l a_l^2 \gg \rho_g a_g^2$ 인 한계를 넣어보자. 분자는 $-\rho_l a_l^2$ , 분모는 $\rho_l a_l^2(1-\psi)$ 에 수렴하므로

K \to -\psi, \qquad \psi + K \to 0

이다. 물리적으로는 당연하다. 물은 거의 압축되지 않으니 물의 체적은 이송만 되어야 한다. $(\psi+K)\nabla\cdot u$ 를 0으로 만드는 것이 바로 그 거동이다.

그렇다면 두 항을 아예 쓰지 않으면 되지 않을까. 순수 액체 셀은 이유가 다르다. $\psi=1$ 이고 $K=0$ 이라 $(\psi+K)\nabla\cdot u = \nabla\cdot u$ 다. 액체도 자기 압축률로는 눌려야 하고, 그 항이 그 일을 한다. 상쇄가 완전한지도 체적탄성계수 비에 달려 있다. 비가 $10$ 근처면 $\psi=0.5$ 에서 $\psi+K=0.0909$ 다. $0$ 이 아니다.

물 기둥을 10 bar로 눌러보면 얼마나 어긋나나#

오일러 격자에서 직접 재자. 구간 $[0,1]$ 의 왼쪽 절반을 물, 오른쪽 절반을 공기로 채운다. 왼쪽은 벽( $u=0$ ), 오른쪽은 공기가 들어오는 열린 경계다. 압력은 공간적으로 균일하게 1 bar에서 10 bar로 선형 상승시킨다. Stiffened-gas 상태방정식( $\gamma_l=4.4$ , $\Pi_l=6\times10^8$ Pa)을 쓰고, 면 속도는 위에서 유도한 혼합물 압축률로 적분한다.

물은 오른쪽 경계에 닿지 않으므로 액체 질량이 보존된다. 따라서 정확한 액체 체적은 $\psi_0\,\rho_l(p_0)/\rho_l(p)$ 로 바로 나온다. 비교할 기준이 공짜로 생긴다.

import math
 
P0, PEND = 1.0e5, 1.0e6                    # 1 bar -> 10 bar
GAM_L, PI_L, RHO_L0 = 4.4, 6.0e8, 1000.0   # water  (stiffened gas)
GAM_G, PI_G, RHO_G0 = 1.4, 0.0,   1.2      # air
 
def stiffened_rho(p, gam, pi, rho0):
    return rho0 * ((p + pi) / (P0 + pi)) ** (1.0 / gam)
 
def stiffened_bulk(p, gam, pi):
    return gam * (p + pi)                  # rho*a^2
 
def wood_speed(psi, p):
    rho_m = psi * stiffened_rho(p, GAM_L, PI_L, RHO_L0) \
          + (1 - psi) * stiffened_rho(p, GAM_G, PI_G, RHO_G0)
    comp  = psi / stiffened_bulk(p, GAM_L, PI_L) \
          + (1 - psi) / stiffened_bulk(p, GAM_G, PI_G)
    return math.sqrt(1.0 / (rho_m * comp)), rho_m
 
def compressibility_factor(psi, p, mode):
    bl = stiffened_bulk(p, GAM_L, PI_L)
    bg = stiffened_bulk(p, GAM_G, PI_G)
    k = (bg - bl) * psi * (1 - psi) / (bg * psi + bl * (1 - psi))
    return 0.0 if mode == 'zero' else (-k if mode == 'flip' else k)
 
def squeeze_column(mode, ncell, nstep, psi0=0.5):
    """왼쪽 벽 u=0, 오른쪽은 순수 기체 유입, p(t)는 1->10 bar 선형"""
    dx, dt = 1.0 / ncell, 1.0 / nstep
    psi = [1.0 if (i + 0.5) * dx < psi0 else 0.0 for i in range(ncell)]
    for n in range(nstep):
        p    = P0 + (PEND - P0) * (n * dt)
        pdot = PEND - P0
        bl   = stiffened_bulk(p, GAM_L, PI_L)
        bg   = stiffened_bulk(p, GAM_G, PI_G)
        u = [0.0] * (ncell + 1)            # 혼합물 압축률 -> du/dx
        for i in range(ncell):
            u[i + 1] = u[i] - dx * pdot * (psi[i] / bl + (1 - psi[i]) / bg)
        nxt = psi[:]
        for i in range(ncell):
            ur, ul = u[i + 1], u[i]
            right  = psi[i + 1] if i + 1 < ncell else 0.0
            fr = ur * (psi[i] if ur > 0 else right)
            fl = (ul * (psi[i - 1] if ul > 0 else psi[i])) if i > 0 else 0.0
            src = (psi[i] + compressibility_factor(psi[i], p, mode)) * (ur - ul) / dx
            nxt[i] = min(1.0, max(0.0, psi[i] - dt * (fr - fl) / dx + dt * src))
        psi = nxt
    exact = psi0 * stiffened_rho(P0, GAM_L, PI_L, RHO_L0) \
                 / stiffened_rho(PEND, GAM_L, PI_L, RHO_L0)
    return sum(psi) * dx, exact
 
print("[1] stiffened-gas properties")
for p in (P0, PEND):
    rl = stiffened_rho(p, GAM_L, PI_L, RHO_L0)
    rg = stiffened_rho(p, GAM_G, PI_G, RHO_G0)
    print(f"  p={p/1e5:4.1f} bar  rho_l={rl:8.3f}  rho_g={rg:6.3f}"
          f"  rho_l*a_l^2={stiffened_bulk(p,GAM_L,PI_L):.3e}"
          f"  rho_g*a_g^2={stiffened_bulk(p,GAM_G,PI_G):.3e}")
 
print("\n[2] K and (psi+K) at 1 bar")
for psi in (0.10, 0.25, 0.50, 0.75, 0.90):
    k = compressibility_factor(psi, P0, 'acid')
    a, _ = wood_speed(psi, P0)
    print(f"  psi={psi:4.2f}  K={k:+9.6f}  psi+K={psi+k:+.3e}  a_wood={a:7.2f} m/s")
 
print("\n[3] liquid-volume change after 1->10 bar compression")
print("  N     exact        ACID         K=0          K flipped    | x exact")
for ncell in (50, 100, 200, 400):
    out = {}
    for mode in ('acid', 'zero', 'flip'):
        v, exact = squeeze_column(mode, ncell, 10 * ncell)
        out[mode] = v - 0.5
    d = exact - 0.5
    print(f"  {ncell:4d}  {d:.4e}  {out['acid']:.4e}  {out['zero']:.4e}  "
          f"{out['flip']:.4e} |  {out['zero']/d:5.3f}  {out['flip']/d:5.3f}")

[1] stiffened-gas properties
  p= 1.0 bar  rho_l=1000.000  rho_g= 1.200  rho_l*a_l^2=2.640e+09  rho_g*a_g^2=1.400e+05
  p=10.0 bar  rho_l=1000.341  rho_g= 6.215  rho_l*a_l^2=2.644e+09  rho_g*a_g^2=1.400e+06
 
[2] K and (psi+K) at 1 bar
  psi=0.10  K=-0.099994  psi+K=+5.891e-06  a_wood=  39.23 m/s
  psi=0.25  K=-0.249982  psi+K=+1.767e-05  a_wood=  27.28 m/s
  psi=0.50  K=-0.499947  psi+K=+5.302e-05  a_wood=  23.65 m/s
  psi=0.75  K=-0.749841  psi+K=+1.590e-04  a_wood=  27.32 m/s
  psi=0.90  K=-0.899523  psi+K=+4.770e-04  a_wood=  39.43 m/s
 
[3] liquid-volume change after 1->10 bar compression
  N     exact        ACID         K=0          K flipped    | x exact
    50  -1.7027e-04  -1.7027e-04  -3.1695e-04  -7.1967e-04 |  1.861  4.227
   100  -1.7027e-04  -1.7027e-04  -3.1559e-04  -7.0606e-04 |  1.853  4.147
   200  -1.7027e-04  -1.7027e-04  -3.1256e-04  -6.7737e-04 |  1.836  3.978
   400  -1.7027e-04  -1.7027e-04  -3.0670e-04  -6.2727e-04 |  1.801  3.684

K = 0, 부호 반전, 체적탄성계수 비 슬라이더가 액체 기둥이 줄어드는 속도를 바꾼다.

ψ₀ (initial liquid fraction) = 0.50final p/p₀ = 10log₁₀(ρ_l a_l² / ρ_g a_g²) = 4.3

p/p₀ = 1.00  |  ψ = 0.50000  |  K = 0.00000  |  ψ+K = 5.00e-1
V_liquid: scheme 0.50000 vs exact 0.50000  →  +0.00%

K = 0을 켜고 재생하면 파란 블록이 정확해보다 빨리 줄어든다. 비를 $1.0$ 까지 내리면 세 기둥이 겹치고, $\psi+K$ 표시가 $0$ 에서 멀어진다.

압력을 10배 올리는 동안 물이 실제로 줄어든 체적은 $1.70\times10^{-4}$ 다. 초기 체적의 0.034%다. ACID 종결은 이 값을 소수점 아래 다섯 자리까지 맞춘다. $K=0$ 은 $3.13\times10^{-4}$ 를 깎았다. 물리적 압축량의 1.84배다.

0.034%와 0.062%의 차이는 작아 보인다. 그러나 비교 대상은 전체 체적이 아니라 물리적 압축량이다. 그 양의 84%를 수치 오차가 덧붙인 것이다. 기포 붕괴나 캐비테이션처럼 체적 변화 자체가 답인 문제에서는 이 비율이 답을 가른다.

표 [2]는 대칭에 가깝다. $\psi=0.10$ 에서 $K=-0.099994$ , $\psi+K=5.891\times10^{-6}$ , 음속 $39.23$ m/s다. $\psi=0.90$ 에서는 $K=-0.899523$ , $\psi+K=4.770\times10^{-4}$ , 음속 $39.43$ m/s다. 양 끝 음속 차는 $0.20$ m/s인데, $\psi+K$ 는 약 $81$ 배 벌어진다. 바닥은 $\psi=0.50$ 의 $23.65$ m/s다. $39.23$ 보다 $15.58$ m/s, $40\%$ 가까이 느리다. 체적탄성계수 비 $2.64\times10^{9}/1.40\times10^{5}$ 는 $18857$ 이다. 그 비가 $\psi=0.5$ 의 $K$ 를 $-0.5$ 바로 옆, 표의 $-0.499947$ 까지 붙인다. $\psi+K$ 는 $+5.302\times10^{-5}$ 다. 표 [3]에서 $K=0$ 배율은 $N=50$ 의 $1.861$ 에서 $N=400$ 의 $1.801$ 로 $0.060$ 줄어든다. 부호를 뒤집은 배율은 $4.227$ 에서 $3.684$ 로 $0.543$ 줄어든다. 같은 $8$ 배 조임인데 부호 쪽 감소가 $0.060$ 의 약 $9$ 배다. $N=400$ 에서도 $1.801$ 과 $3.684$ 는 $1$ 에서 멀다. $N=50$ 에서 $3.1695\times10^{-4}$ 와 $1.7027\times10^{-4}$ 의 차는 $1.4668\times10^{-4}$ 이고, 본문의 $1.467\times10^{-4}$ 로 반올림된다. $N=400$ 은 $3.0670\times10^{-4}$ 에서 $1.7027\times10^{-4}$ 를 빼 $1.3643\times10^{-4}$ 다. $10$ bar에서 액체 밀도는 $1000.000$ 에서 $1000.341$ 로 $0.341$ 오르고, 기체는 $1.200$ 에서 $6.215$ 로 뛴다. 없는 $K$ 는 그 작은 액체 압축의 $1.86$ 배를 계면에 청구한다.

격자를 8배 조이면 그 오차는 사라지나#

표 [3]의 마지막 두 열이 그 답이다. $N$ 을 50에서 400으로 8배 늘리는 동안 $K=0$ 의 비는 1.861에서 1.801로 움직였다. 오차 크기로는 $1.467\times10^{-4}$ 에서 $1.364\times10^{-4}$ , 7% 감소다. 수렴이라고 부를 수 없다.

이유는 $K=0$ 이 이산화 오차가 아니기 때문이다. $K=0$ 으로 두면 푸는 방정식 자체가 $\mathrm{D}\psi/\mathrm{D}t=0$ 으로 바뀐다. 다른 방정식의 해로 수렴하는 것이다. $\Delta x\to0$ 에서 사라질 이유가 없다.

남은 7%는 계면이 번진 띠의 모양이 격자에 따라 조금 달라지기 때문이다. 오차의 크기는 그 띠 안에서 $\psi\,\nabla\cdot u$ 를 적분한 양이고, 1차 상류차분의 수치 확산은 격자를 조여도 띠를 셀 개수로 넓힌다. 띠의 물리적 두께가 거의 그대로이므로 적분값도 거의 그대로다. 그 띠를 한 스텝에 얼마나 밀 수 있는지는 VOF 계면 이류의 CFL 상한이 따로 정한다.

이 성질은 나쁜 소식이자 좋은 소식이다. 격자 수렴성 시험으로는 이 버그를 못 찾는다. 반대로, 격자 두 개만 돌려 비가 일정하게 남으면 종결이 틀렸다는 신호다.

부호를 반대로 적으면 어떤 모양으로 틀리나#

세 번째 열이 부호를 뒤집은 결과다. $K\to-K$ 이면 $\psi+K \to 2\psi$ 가 되고, 압축량은 3.68~4.23배로 뛴다. 방향도 뒤집힌다. 압축하는데 액체 체적률이 내려가므로, 화면에서는 기포가 자라는 것처럼 보인다.

문제는 이쪽도 그럴듯한 그림을 준다는 점이다. 계면은 여전히 매끈하고, 체적률은 여전히 $[0,1]$ 안에 있고, 발산하지도 않는다. 단상 한계만 확인해도 못 잡는다. $\psi\to0$ 과 $\psi\to1$ 에서 $K$ 는 두 부호 모두 0이기 때문이다.

그래서 논문 식을 코드에 옮길 때는 $\psi$ 가 어느 상의 체적률인지 주석으로 못 박고, 위의 1차원 압축 테스트 하나로 부호를 고정하는 편이 빠르다. 테스트는 격자 50개로 1초면 돌아간다.

같은 가정이 엔탈피에서 한 번 더 쓰이는 자리#

논문의 ACID(Acoustically-Conservative Interface Discretisation)는 이 등압 종결을 밀도와 엔탈피에 각각 적용한다. 밀도는 쉽다. 부분밀도를 선형 보간하면

\rho = \rho_g + \psi(\rho_l - \rho_g)

이고, 이것이 등압 종결과 수학적으로 동일하다. 질량과 운동량 보존이 그대로 따라온다.

엔탈피는 다르다. $h$ 는 선형계에서 직접 푸는 1차 변수이므로 밀도처럼 보간식을 갈아끼울 수 없다. 그래서 ACID 가정을 적용한 타겟 면 엔탈피를 따로 계산한다.

\hat{h}_f = \frac{1}{\tilde{\rho}_f}\left[\rho_U h_U + \xi_f L_f (\rho_D h_D - \rho_U h_U)\right]

$U$ 와 $D$ 는 상류·하류 셀, $L_f$ 는 보간 가중치, $\xi_f$ 는 제한자다. 그리고 현재 값과의 차이를 지연 보정항으로 에너지 방정식 우변에 더한다.

\delta h_f = \hat{h}_f - \tilde{h}_f

$K$ 와 $\delta h_f$ 는 같은 가정의 두 얼굴이다. 체적률 식에서는 대수적으로 닫혀 계수 하나로 떨어지고, 1차 변수인 엔탈피에서는 닫히지 않아 명시적 소스로 미뤄진다. 미뤄진 쪽은 반복 내에서 수렴해야 하므로 비선형 반복 횟수에 영향을 준다. 등압으로 한 압력에 붙이는 마무리는 두 압력을 풀고 한 압력으로 이완하는 경로의 끝과 같다. 차이는 이쪽이 그 이완을 $K$ 하나로 접는다는 점이다.

Denner가 리만 솔버를 버릴 수 있었던 이유#

밀도 기반 해법은 계면을 리만 문제로 풀어 접촉 불연속면을 유지한다. 압력기반 완전연성 해법은 $\rho$ 를 해석 변수에서 빼고 $u$ , $p$ , $h$ 를 한 선형계로 묶는다. 동일 위치 격자에서 생기는 압력-속도 분리는 MWI로 막는다. 극심한 밀도비에서는 조화 평균 면 밀도를 가중치로 쓴다. 접촉면을 재구성으로 지키는 쪽은 더 작은 진폭을 고르는 MUSCL-THINC-BVD다.

리만 솔버를 빼면서 잃는 것은 계면에서의 물리를 자동으로 잡아주던 장치다. 그 자리를 종결 가정으로 직접 메워야 한다. $K$ 와 $\delta h_f$ 가 그 메운 부분이다.

혼합물 음속 23.65 m/s를 음해적으로 건너뛰면서 접촉 불연속면이 번지지 않게 하려면 그 종결이 정확해야 한다. 음향을 음해로 넘기면 CFL 청구서가 수송 쪽으로 옮겨 간다. 음향과 수송을 나눠 CFL을 푸는 분리가 그 자리다. 1.86배는 그 종결이 하는 일의 크기를 재어본 숫자다.