아무 일도 일어나면 안 되는 문제에서 압력이 21% 튀었다

비보존항의 이산화는 자유롭게 고를 수 있는 것이 아니다. 보존 플럭스를 먼저 정하면 균일 유동 조건이 나머지 항의 계수를 하나로 못박는다.

아무 일도 일어나면 안 되는 문제가 첫 번째 테스트였다#

다성분 압축성 솔버의 첫 검증 케이스를 골랐다. 두 기체가 접촉면 하나를 사이에 두고 맞닿아 있다. 압력은 어디서나 100 kPa, 속도는 어디서나 100 m/s. 두 기체의 밀도와 비열비만 다르다.

정확해는 지루하다. 접촉면이 오른쪽으로 흘러갈 뿐, 압력과 속도는 끝까지 균일하다. 격자가 거칠어도, 시간간격을 크게 잡아도 마찬가지다. 이 문제에는 풀어야 할 파동이 없다.

20스텝 뒤 계면 셀의 압력이 78.6 kPa로 내려갔다. 균일값에서 21% 벗어난 값이다. 속도는 정확히 100 m/s를 유지했고, 질량과 에너지도 기계 정밀도까지 보존됐다. 플럭스는 멀쩡했다. 틀린 것은 플럭스 옆에 서 있던 항이었다.

아래 시뮬레이션에서 그 진동을 직접 만들어보자.

gamma_21.667rho_1/rho_28.0CFL0.50

step 0A  max|P-P0| 0.00 kPa (0.0%)B  max|P-P0| 0.0e+0 Pa

Drag gamma_2 toward 1.400: the red curve flattens onto the green one, because with equal gamma the mixture rule is no longer nonlinear. Push it away, or raise the density ratio, and the red curve dips at the two interfaces while the green line never leaves 100 kPa. CFL changes how fast the damage accumulates, not whether it appears.

두 곡선은 같은 보존변수를 같은 풍상 연산자로 옮긴다. 다른 것은 열역학 변수를 태우는 방식뿐이다. gamma_2를 1.400 쪽으로 끌면 붉은 곡선이 초록 곡선 위로 내려앉는다. 반대로 밀거나 밀도비를 올리면 계면 두 곳에서 붉은 곡선이 아래로 꺼진다. 압력이 균일하게 유지되는지는 스킴의 정확도가 아니라 이 선택 하나로 갈린다.

이 계산에서 정확히 무엇이 어긋났는가#

압력은 보존변수가 아니다. 상태방정식으로 되찾는 값이다. 비열비 $\gamma$ 가 공간에 따라 변하면

P = \frac{\rho E - \tfrac{1}{2}\rho u^{2}}{\Gamma}, \qquad \Gamma \equiv \frac{1}{\gamma - 1}

여기서 $\rho E$ 는 단위 부피당 총에너지, $\Gamma$ 는 비열비를 담은 열역학 변수다.

균일 압력·속도에서 총에너지는 $\rho E = P\Gamma + \tfrac{1}{2}\rho u^{2}$ 이다. 풍상 이류는 선형이므로 $\rho E$ 의 갱신값은 $P \cdot (\Gamma$ 의 풍상 이류 $)$ 에 운동에너지 항을 더한 값이 된다. 그러면 새 압력은 이렇게 정리된다.

P^{n+1}_j = P^{n} \, \frac{\mathcal{U}(\Gamma)_j}{\Gamma^{n+1}_j}

$\mathcal{U}$ 는 $\rho E$ 에 실제로 적용된 그 풍상 연산자다.

조건은 분수 하나로 남는다. $\Gamma^{n+1}$ 을 $\mathcal{U}(\Gamma)$ 와 같은 값으로 만들어야 압력이 제자리에 머문다. 즉 $\Gamma$ 를 옮기는 방식이 에너지를 옮기는 방식과 대수적으로 맞물려야 한다.

내가 처음 쓴 코드는 질량분율 $Y$ 를 보존형으로 옮기고 혼합 규칙으로 $\Gamma$ 를 되돌렸다. $\Gamma(Y)$ 는 $Y$ 에 대해 선형이지만, $Y^{n+1} = \mathcal{U}(\rho Y)/\mathcal{U}(\rho)$ 는 밀도가 변하는 곳에서 $\mathcal{U}(Y)$ 가 아니다. 분수 두 개가 어긋나면 그 차이가 곧 압력 오차다.

Abgrall이 1996년에 세운 요구사항 한 줄#

Re와 Abgrall이 약압축성 다성분 모델을 세울 때 인용한 기준은 이렇게 한 문장이다. "압력과 속도가 균일한 2상 유동은 시간이 지나도 같은 변수에 대해 균일하게 유지되어야 한다." 논문은 이것을 압력 비교란 조건(pressure non-disturbance condition), 또는 Abgrall의 기준이라 부른다.

이 문장이 특이한 이유는 정확도에 대한 요구가 아니라는 점이다. 1차든 5차든 상관없다. 안정성 조건도 아니다. CFL을 낮춰도 진동은 그대로 남는다. 이것은 이산화들 사이의 대수적 정합에 대한 요구다. 한 방정식에 어떤 연산자를 썼는지가 다른 방정식에서 허용되는 연산자를 결정한다.

같은 사정은 두 유체 모델의 쌍곡성과 계면 압력에서 계면 압력 항을 고를 때도 나타난다. 거기서는 고윳값이 실수가 되도록 항을 골랐고, 여기서는 균일 유동이 유지되도록 항을 고른다. 두 경우 모두 "물리적으로 그럴듯한 이산화"가 여러 개 있고 그중 하나만 살아남는다.

Python으로 같은 격자에 두 이산화를 나란히 세웠다#

100셀 주기 격자에 접촉면 하나를 놓고, 두 방식을 같은 스텝 수만큼 돌렸다. 보존변수 $\rho$ , $\rho u$ , $\rho E$ 는 두 방식 모두 동일한 1차 풍상으로 갱신한다. 차이는 열역학 변수 하나뿐이다.

G1, G2 = 1.4, 1.667      # 두 기체의 비열비
P0, U0 = 1.0e5, 100.0    # 균일 압력 [Pa], 균일 속도 [m/s]
R1, R2 = 1.0, 0.125      # 두 기체의 밀도 [kg/m^3]
 
 
def gamma_var(y):
    """질량분율 y -> 1/(gamma-1). y에 대해 선형인 혼합 규칙."""
    return y / (G1 - 1.0) + (1.0 - y) / (G2 - 1.0)
 
 
def advect_upwind(q, lam):
    """u>0 1차 풍상 이류. 좌측 유입 셀은 고정."""
    return [q[0]] + [q[j] - lam * (q[j] - q[j - 1]) for j in range(1, len(q))]
 
 
def initial_state(n):
    x = [(j + 0.5) / n for j in range(n)]
    y = [1.0 if xi < 0.3 else 0.0 for xi in x]
    rho = [R1 if xi < 0.3 else R2 for xi in x]
    return x, y, rho
 
 
def step_massfraction_closure(rho, mom, ene, ry, lam):
    """rho*Y를 보존형으로 옮기고 혼합 규칙으로 gamma를 되돌린다."""
    rho_n = advect_upwind(rho, lam)
    mom_n = advect_upwind(mom, lam)
    ene_n = advect_upwind(ene, lam)
    ry_n = advect_upwind(ry, lam)
    y_n = [ry_n[j] / rho_n[j] for j in range(len(rho_n))]
    p_n = [(ene_n[j] - 0.5 * mom_n[j] ** 2 / rho_n[j]) / gamma_var(y_n[j])
           for j in range(len(rho_n))]
    return rho_n, mom_n, ene_n, ry_n, p_n
 
 
def step_gammavar_transport(rho, mom, ene, gv, lam):
    """1/(gamma-1)을 비보존(이류) 형태로 같은 풍상 연산자에 태운다."""
    rho_n = advect_upwind(rho, lam)
    mom_n = advect_upwind(mom, lam)
    ene_n = advect_upwind(ene, lam)
    gv_n = advect_upwind(gv, lam)
    p_n = [(ene_n[j] - 0.5 * mom_n[j] ** 2 / rho_n[j]) / gv_n[j]
           for j in range(len(rho_n))]
    return rho_n, mom_n, ene_n, gv_n, p_n
 
 
def run_interface_advection(n=100, steps=60, cfl=0.5):
    x, y0, rho0 = initial_state(n)
    gv0 = [gamma_var(v) for v in y0]
 
    rho_a = list(rho0)
    mom_a = [r * U0 for r in rho0]
    ene_a = [P0 * gv0[j] + 0.5 * rho0[j] * U0 ** 2 for j in range(n)]
    ry_a = [rho0[j] * y0[j] for j in range(n)]
 
    rho_b, mom_b, ene_b = list(rho_a), list(mom_a), list(ene_a)
    gv_b = list(gv0)
 
    hist, p_a, p_b = [], None, None
    for k in range(1, steps + 1):
        rho_a, mom_a, ene_a, ry_a, p_a = step_massfraction_closure(
            rho_a, mom_a, ene_a, ry_a, cfl)
        rho_b, mom_b, ene_b, gv_b, p_b = step_gammavar_transport(
            rho_b, mom_b, ene_b, gv_b, cfl)
        if k % 20 == 0:
            ea = max(abs(v - P0) for v in p_a)
            eb = max(abs(v - P0) for v in p_b)
            eu = max(abs(mom_a[j] / rho_a[j] - U0) for j in range(n))
            hist.append((k, ea, eb, eu))
    return x, p_a, p_b, hist
 
 
if __name__ == "__main__":
    x, p_a, p_b, hist = run_interface_advection()
    print("step |  max|P-P0| mixrule |  max|P-P0| Gamma-adv |  max|u-U0| mixrule")
    for k, ea, eb, eu in hist:
        print("%4d | %16.2f | %19.2e | %16.3e" % (k, ea, eb, eu))
 
    j = max(range(len(p_a)), key=lambda i: abs(p_a[i] - P0))
    print("\nworst cell x=%.3f  P=%.1f Pa  (uniform value %.0f Pa)" % (x[j], p_a[j], P0))
    print("relative error:  mixrule %.2f%%   Gamma-adv %.1e%%"
          % (max(abs(v - P0) for v in p_a) / P0 * 100,
             max(abs(v - P0) for v in p_b) / P0 * 100))

step |  max|P-P0| mixrule |  max|P-P0| Gamma-adv |  max|u-U0| mixrule
  20 |         21433.25 |            2.91e-11 |        0.000e+00
  40 |         21587.18 |            4.37e-11 |        0.000e+00
  60 |         21442.28 |            4.37e-11 |        2.842e-14
 
worst cell x=0.635  P=78557.7 Pa  (uniform value 100000 Pa)
relative error:  mixrule 21.44%   Gamma-adv 4.4e-14%

세 줄만 읽으면 된다. 혼합 규칙 쪽은 21.4 kPa에서 시작해 스텝을 늘려도 줄지 않는다. $\Gamma$ 이류 쪽은 $10^{-11}$ Pa, 상대오차로 $10^{-14}$ %다. 배정밀도의 바닥이다.

속도 열이 0인 것도 중요하다. 운동량 방정식은 처음부터 끝까지 옳게 풀렸다. 격자를 두 배로 조밀하게 해도 21%는 21%로 남는다. 수렴하지 않는 오차이기 때문이다.

논문이 $H_u$ 를 유도한 방식 — 스킴을 먼저 정하고 항을 나중에 맞춘다#

Re와 Abgrall의 Baer–Nunziato형(두 상이 각자 속도·압력을 갖는 7방정식 계열) 모델에는 부피분율 방정식이 따로 있다. 이 방정식은 보존형이 아니다.

\frac{\partial \alpha_i}{\partial t} + u_I \frac{\partial \alpha_i}{\partial x} = 0

$\alpha_i$ 는 상 $i$ 의 부피분율, $u_I$ 는 계면 속도다. 논문은 이 항을 $H_u(\alpha_i, u_I)$ 라는 이름의 이산 연산자로 두고, 그 형태를 가정하지 않는다. 대신 조건에서 끌어낸다.

질량 방정식 $(\alpha_i\rho_i)$ 는 Rusanov 플럭스로 이미 정해져 있다. 밀도와 속도가 균일한 상태를 넣으면 그 플럭스는 $\rho_i \times (\alpha_i$ 의 Rusanov 플럭스 $)$ 로 인수분해된다. 갱신 후에도 $\rho_i = (\alpha_i\rho_i)^{n+1}/\alpha_i^{n+1}$ 가 그대로 남으려면, 분모 $\alpha_i^{n+1}$ 이 분자와 똑같은 플럭스 차분으로 갱신돼야 한다. 그래서 $H_u$ 는 선택이 아니라 결과다.

H_u(\alpha_i, u_I)_j = \frac{1}{2}\Big[ \big(\alpha_{j+1} - \alpha_{j-1}\big) u_{I,j} - \big|u_{I,j}\big| \big(\alpha_{j+1} - 2\alpha_j + \alpha_{j-1}\big) \Big]

앞의 항은 중심차분, 뒤의 항은 $|u_I|$ 가 붙은 확산이다. 둘의 합은 정확히 $\alpha$ 에 대한 Rusanov 플럭스 차분이다. $u_I > 0$ 이면 $\alpha_j - \alpha_{j-1}$ 만 남는다. 순수 풍상이다.

논문은 압력 방정식에서도 같은 $H_u$ 를 재사용한다. $u^{*}\partial P/\partial x$ 를 $\partial(Pu^{*})/\partial x - P\partial u^{*}/\partial x$ 로 쪼갠 뒤, 남은 비보존항을 질량 방정식과 같은 연산자에 태운다. 방정식마다 다른 이산화를 쓰면 방금 맞춰둔 정합이 깨지기 때문이다.

세 점 스텐실의 계수를 손으로 돌려본다#

$H_u$ 의 두 항 중 확산 항에 가중치 $\theta$ 를 붙이면 한 다이얼로 두 극단을 오갈 수 있다. $\theta = 1$ 이면 논문의 $H_u$ , $\theta = 0$ 이면 $u_I \partial\alpha/\partial x$ 의 순수 중심차분이다.

theta1.00CFL0.40

step 0weights [-1.00, 1.00, 0.00]max|rho - 850| 0.0e+0 kg/m^3

Slide theta down from 1 and watch the downwind weight come back to life: the moment the third cell rejoins the stencil, the recovered density leaves 850 and never returns. Flip u_I and theta = 1 still holds, because the flux difference follows the sign of the interface velocity while the centred stencil does not.

theta를 1에서 내리면 하류 셀의 가중치가 0에서 되살아난다. 그 순간 되찾은 밀도가 850 kg/m³을 떠나 돌아오지 않는다. u_I 부호를 뒤집어도 $\theta = 1$ 은 버틴다. 플럭스 차분은 계면 속도의 부호를 따라가지만 중심 스텐실은 그러지 않기 때문이다.

숫자로도 확인했다. 균일 밀도 850 kg/m³, 계면 하나, 같은 격자에서 $\alpha$ 만 두 연산자로 갱신했다.

RHO, UI, N, LAM = 850.0, 1.0, 80, 0.4   # 균일 밀도 [kg/m^3], 계면속도, 셀 수, u*dt/dx
 
 
def alpha_profile():
    """계면을 사이에 둔 부피분율. 0.02 <-> 0.98 사이를 세 셀에 걸쳐 잇는다."""
    a = []
    for j in range(N):
        if j < 30:
            a.append(0.98)
        elif j < 33:
            a.append(0.98 - 0.32 * (j - 29))
        else:
            a.append(0.02)
    return a
 
 
def rusanov_flux(q, j, vel):
    """셀 j와 j+1 사이의 Rusanov 수치 플럭스 (주기 경계)."""
    ql, qr = q[j % N], q[(j + 1) % N]
    return 0.5 * (qr + ql) * vel - 0.5 * abs(vel) * (qr - ql)
 
 
def hu_upwind(a, j):
    """논문 식 (10)의 비보존 연산자: alpha에 대한 Rusanov 플럭스 차분."""
    return rusanov_flux(a, j, UI) - rusanov_flux(a, j - 1, UI)
 
 
def hu_central(a, j):
    """중심차분으로 u_I * d(alpha)/dx 를 그대로 이산화한 버전."""
    return 0.5 * UI * (a[(j + 1) % N] - a[(j - 1) % N])
 
 
def march(op, steps):
    """alpha*rho는 Rusanov, alpha는 op로 전진시키고 되찾은 rho를 본다."""
    a = alpha_profile()
    ar = [RHO * v for v in a]
    for _ in range(steps):
        ar = [ar[j] - LAM * (rusanov_flux(ar, j, UI) - rusanov_flux(ar, j - 1, UI))
              for j in range(N)]
        a = [a[j] - LAM * op(a, j) for j in range(N)]
    return max(abs(ar[j] / a[j] - RHO) for j in range(N))
 
 
if __name__ == "__main__":
    print("steps |  upwind H_u [kg/m^3] |  centred [kg/m^3]")
    for s in (10, 40, 120):
        print("%5d | %20.2e | %17.4f" % (s, march(hu_upwind, s), march(hu_central, s)))
 
    a = alpha_profile()
    lhs = hu_upwind(a, 31)
    rhs = 0.5 * ((a[32] - a[30]) * UI - abs(UI) * (a[32] - 2 * a[31] + a[30]))
    print("\ncell 31:  flux difference %.6f   paper Eq.(10) %.6f   gap %.1e"
          % (lhs, rhs, abs(lhs - rhs)))

steps |  upwind H_u [kg/m^3] |  centred [kg/m^3]
   10 |             2.27e-13 |         2553.2027
   40 |             4.55e-13 |         6697.6224
  120 |             5.68e-13 |         1206.6015
 
cell 31:  flux difference -0.320000   paper Eq.(10) -0.320000   gap 0.0e+00

풍상 $H_u$ 는 120스텝을 돌아도 $10^{-13}$ kg/m³에 머문다. 중심차분은 10스텝 만에 2,553 kg/m³을 벗어난다. 120스텝에서 값이 다시 1,207로 내려온 것은 회복이 아니라 발산이다. $\alpha$ 가 0 근처를 오가면서 나눗셈이 아무 값이나 뱉는 구간에 들어갔다는 뜻이다.

마지막 줄은 유도가 맞았다는 확인이다. $\alpha$ 에 대한 Rusanov 플럭스 차분과 논문 식 (10)의 닫힌 형태가 셀 31에서 자릿수까지 같다. 두 표현은 대수적으로 같은 것이다.

이 조건이 정확도와 무관하다는 말의 의미#

한 가지를 분명히 해두자. $H_u$ 를 쓴다고 해가 정확해지는 것은 아니다. 1차 풍상은 여전히 계면을 뭉갠다. 위 그림에서도 $\alpha$ 의 계단은 스텝마다 두꺼워진다.

Abgrall 기준이 보장하는 것은 다른 종류다. 틀리더라도 물리적으로 말이 되는 방향으로 틀린다. 계면이 번지는 것은 수치 확산이고, 그것은 격자를 조밀하게 하면 줄어든다. 균일 압력장에 21 kPa짜리 스파이크가 서는 것은 수치 확산이 아니다. 그것은 격자를 조밀하게 해도 줄지 않고, 상태방정식이 강성(stiff)일수록 커지며, 종종 음의 압력으로 계산을 멈춘다.

암시적 표면장력이 시간간격을 여는 폭을 다룰 때도 같은 구분이 나왔다. 어떤 제약은 정확도를 사면 풀리고, 어떤 제약은 사도 풀리지 않는다. Abgrall 기준은 후자다. 고차 스킴으로 올라가도 이 조건은 따로 다시 맞춰야 한다.

다음에 계면에서 압력이 튀면 어디부터 보는가#

내가 이 문제를 다시 만나면 순서는 이렇다.

먼저 균일 압력·속도 문제를 돌린다. 파동이 없는 문제에서 압력 변동이 0이 아니면 플럭스를 볼 필요가 없다. 답은 비보존항이나 상태방정식 복원 쪽에 있다.

그다음 격자를 두 배로 조밀하게 해서 같은 값을 잰다. 진동이 절반으로 줄면 수치 확산 문제다. 그대로면 정합성 문제다. 이 한 번의 실행이 두 원인을 갈라준다.

마지막으로 방정식마다 비보존항에 어떤 스텐실을 썼는지 나란히 적어본다. 질량은 Rusanov, 부피분율은 중심차분, 에너지는 또 다른 무엇이었다면 거기가 원인이다. 논문이 $H_u$ 하나를 세 방정식에 돌려쓴 것은 코드를 줄이려는 취향이 아니었다.

참고 문헌

B. Re, R. Abgrall, Non-equilibrium Model for Weakly Compressible Multi-component Flows: the Hyperbolic Operator, arXiv:1911.00270 — §2.2 The discretization
R. Abgrall, How to Prevent Pressure Oscillations in Multicomponent Flow Calculations: A Quasi Conservative Approach, J. Comput. Phys. 125 (1996)

아무 일도 일어나면 안 되는 문제에서 압력이 21% 튀었다 — Abgrall 기준과 비보존항

아무 일도 일어나면 안 되는 문제가 첫 번째 테스트였다#

이 계산에서 정확히 무엇이 어긋났는가#

Abgrall이 1996년에 세운 요구사항 한 줄#

Python으로 같은 격자에 두 이산화를 나란히 세웠다#

논문이 $H_u$ 를 유도한 방식 — 스킴을 먼저 정하고 항을 나중에 맞춘다#

세 점 스텐실의 계수를 손으로 돌려본다#

이 조건이 정확도와 무관하다는 말의 의미#

다음에 계면에서 압력이 튀면 어디부터 보는가#

관련

아무 일도 일어나면 안 되는 문제에서 압력이 21% 튀었다 — Abgrall 기준과 비보존항

아무 일도 일어나면 안 되는 문제가 첫 번째 테스트였다#

이 계산에서 정확히 무엇이 어긋났는가#

Abgrall이 1996년에 세운 요구사항 한 줄#

Python으로 같은 격자에 두 이산화를 나란히 세웠다#

논문이 HuH_uHu​를 유도한 방식 — 스킴을 먼저 정하고 항을 나중에 맞춘다#

세 점 스텐실의 계수를 손으로 돌려본다#

이 조건이 정확도와 무관하다는 말의 의미#

다음에 계면에서 압력이 튀면 어디부터 보는가#

관련

논문이 $H_u$ 를 유도한 방식 — 스킴을 먼저 정하고 항을 나중에 맞춘다#