격자 Péclet를 2로 올리자 해가 -0.33까지 내려갔다 — 최소화할 에너지가 사라진 자리

안정화는 대칭성을 되돌리지 못한다. 대칭 부분을 키워서 이산 최대원리를 대신 사는 것이다.

막대 하나에서는 네 가지 방법이 같은 답을 냈다#

유한요소법 강의노트 1장은 같은 막대 문제를 네 번 푼다. 직접 강성법, 최소 포텐셜 에너지, 가중잔차법, 그리고 갈러킨. 네 번 모두 똑같은 5×5 강성행렬이 나온다. 노트에는 이렇게 적혀 있다. "어떤 것을 사용해도 결과는 별 차이 없이 정확한 값을 구할 수 있다."

그 문장에는 조건이 하나 붙어 있는데, 막대 문제 안에 숨어 있어서 보이지 않는다. 유동 방정식에는 그 조건이 없다. 조건이 깨지면 어떤 보증이 사라지는지, 안정화가 무엇을 되돌리고 무엇을 못 돌리는지를 숫자로 잰다. advection a를 0에서 올리면 분홍 점의 한 바퀴 일이 0.000에서 $2\pi a$ 로 바뀐다.

advection a1.00diffusion s1.00

Set advection a to 0: the pink dot goes all the way around and brings back 0.000 — the field is a gradient, and the amber ball slides straight down the ellipses. Push a up and the last lap returns 6.283, exactly 2πa. Now drag diffusion s across its whole range: the lap total does not budge. No amount of added diffusion buys back a potential.

회색 등에너지선을 가로지르지 않던 노란 공이, a가 커지면 나선을 그리며 거슬러 오른다. diffusion s를 끝까지 밀어도 한 바퀴 일은 그대로다.

최소화할 것이 있으려면 행렬이 대칭이어야 한다#

최소 포텐셜 에너지 정식화는 총 포텐셜 에너지를 절점 변위로 미분해서 0으로 놓는다. 이산화된 형태는 이렇다.

\Pi(\boldsymbol{\phi}) = \frac{1}{2}\,\boldsymbol{\phi}^{T}\mathbf{K}\,\boldsymbol{\phi} - \mathbf{f}^{T}\boldsymbol{\phi}

$\boldsymbol{\phi}$ 는 절점 미지수 벡터, $\mathbf{K}$ 는 강성행렬, $\mathbf{f}$ 는 하중 벡터다. 이제 $i$ 번째 성분으로 미분한다.

\frac{\partial \Pi}{\partial \phi_i} = \frac{1}{2}\sum_{j}\left(K_{ij} + K_{ji}\right)\phi_j - f_i

$K_{ij} + K_{ji}$ 가 나온 것에 주목하자. 이 식이 $\mathbf{K}\boldsymbol{\phi} = \mathbf{f}$ 가 되려면 $K_{ij} = K_{ji}$ 여야 한다. 대칭이 아니면 최소화가 실제로 푸는 것은 대칭화된 행렬 $\tfrac{1}{2}(\mathbf{K}+\mathbf{K}^{T})$ 이고, 그것은 원래 방정식이 아니다.

더 근본적인 이야기도 같은 자리에서 나온다. 잔차장 $\mathbf{r}(\boldsymbol{\phi}) = \mathbf{f} - \mathbf{K}\boldsymbol{\phi}$ 가 어떤 스칼라 함수의 기울기이려면, 그 야코비안이 대칭이어야 한다. 대칭이 아니면 퍼텐셜 함수는 아예 존재하지 않는다. 존재 여부를 확인하는 가장 빠른 방법은 닫힌 경로를 한 바퀴 도는 것이다. 기울기장이라면 한 바퀴 돈 일은 반드시 0이다. 위에서 a를 올렸을 때 공이 등에너지선을 거슬러 오른 것이 그 0이 깨진 자리다.

이류항이 깨뜨리는 양은 정확히 $u$ 다#

1차원 이류-확산 방정식 $u\,\phi' = \epsilon\,\phi''$ 의 약형식을 쓰면 두 항이 성격이 다르다는 것이 드러난다.

a(w, \phi) = \int_0^L \left( \epsilon\,\frac{dw}{dx}\frac{d\phi}{dx} + w\,u\,\frac{d\phi}{dx} \right) dx

$w$ 는 시험함수, $\epsilon$ 은 확산계수, $u$ 는 이류 속도다. 첫 항은 $w$ 와 $\phi$ 를 바꿔도 그대로다 — 대칭이다. 뒤의 항은 부분적분하면 부호가 뒤집힌다. 경계에서 사라지는 시험함수에 대해 $\int w\,u\,\phi'\,dx = -\int \phi\,u\,w'\,dx$ 이므로, 이류항은 순수한 반대칭(skew) 성분이다.

선형 요소로 조립하면 이 구조가 그대로 계수에 남는다. 요소 크기 $h$ , 내부 절점의 행은 이렇다.

K_{i,i-1} = -\frac{\epsilon}{h} - \frac{u}{2}, \qquad K_{i,i} = \frac{2\epsilon}{h}, \qquad K_{i,i+1} = -\frac{\epsilon}{h} + \frac{u}{2}

확산은 양쪽 이웃에 같은 값을 넣고, 이류는 한쪽에 $+u/2$ , 반대쪽에 $-u/2$ 를 넣는다. 그래서 대칭에서 벗어난 양은 정확히 하나의 숫자다.

K_{i,i+1} - K_{i+1,i} = u

격자를 아무리 조밀하게 해도 이 값은 $u$ 다. $h$ 가 들어 있지 않기 때문이다. 유동이 있는 한 최소 포텐셜 에너지 원리는 돌아오지 않는다. 강성행렬이 대칭일 때만 6×6 Voigt 행렬로 접히던 상황과는 사정이 다르다. 이류가 계단 근처에서 2차 스킴을 떨리게 할 때의 처방은 플럭스 제한자의 TVD 이류 쪽이다.

Python으로 같은 격자에 세 가지를 세웠다#

요소 10개, $u = 1$ , 경계조건 $\phi(0)=0$ , $\phi(1)=1$ . 인공확산 계수 $\beta$ 하나로 세 스킴을 만든다. $\epsilon_{\text{eff}} = \epsilon + \beta\,u\,h/2$ 이고, $\beta = 0$ 이면 순수 갈러킨, $\beta = 1$ 이면 완전 상류차분, $\beta = \coth(Pe_h) - 1/Pe_h$ 이면 SUPG다.

import numpy as np
 
def assemble_ad(n, vel, eps, beta):
    """1D advection-diffusion on linear elements. beta = artificial diffusion."""
    h = 1.0 / n
    eps_eff = eps + beta * vel * h / 2.0
    kd = (eps_eff / h) * np.array([[1.0, -1.0], [-1.0, 1.0]])   # diffusion: symmetric
    ka = (vel / 2.0) * np.array([[-1.0, 1.0], [-1.0, 1.0]])     # advection: skew
    K = np.zeros((n + 1, n + 1))
    for e in range(n):
        K[e:e + 2, e:e + 2] += kd + ka
    return K
 
def solve_bvp(K):
    n = K.shape[0] - 1
    A, b = K.copy(), np.zeros(n + 1)
    A[0, :], A[0, 0], b[0] = 0.0, 1.0, 0.0
    A[n, :], A[n, n], b[n] = 0.0, 1.0, 1.0
    return np.linalg.solve(A, b)
 
def exact_ad(x, pe):
    return (np.exp(pe * (x - 1.0)) - np.exp(-pe)) / (1.0 - np.exp(-pe))
 
def skew_ratio(K):
    return np.linalg.norm(K - K.T) / np.linalg.norm(K + K.T)
 
def loop_work(K, m=20000):
    """Work of the residual field -K.phi once around the unit circle in DOF space."""
    t = np.linspace(0.0, 2.0 * np.pi, m, endpoint=False)
    path = np.stack([np.cos(t), np.sin(t)])          # position
    tang = np.stack([-np.sin(t), np.cos(t)])         # dl / dt
    return float(np.sum(np.sum(-(K @ path) * tang, axis=0)) * (2.0 * np.pi / m))
 
n, vel = 10, 1.0
x = np.linspace(0.0, 1.0, n + 1)
print(f"{'Pe_h':>5} {'scheme':>9} {'max|err|%':>10} {'min phi':>9} {'skew/sym':>9} {'loop W':>8}")
for pe_h in [0.5, 1.0, 2.0, 5.0]:
    eps = vel / (2.0 * n * pe_h)
    ex = exact_ad(x, vel / eps)
    for name, beta in [("Galerkin", 0.0), ("upwind", 1.0),
                       ("SUPG", 1.0 / np.tanh(pe_h) - 1.0 / pe_h)]:
        K = assemble_ad(n, vel, eps, beta)
        phi = solve_bvp(K)
        print(f"{pe_h:5.1f} {name:>9} {100 * np.max(np.abs(phi - ex)):10.2f} "
              f"{phi.min():9.4f} {skew_ratio(K):9.4f} {loop_work(K[4:6, 4:6]):8.4f}")
 
K0 = assemble_ad(n, 0.0, 0.1, 0.0)
print(f"\nvel = 0 : skew/sym = {skew_ratio(K0):.2e},  loop W = {loop_work(K0[4:6, 4:6]):.2e}")
print(f"pi * u  = {np.pi * vel:.4f}")

 Pe_h    scheme  max|err|%   min phi  skew/sym   loop W
  0.5  Galerkin       3.45    0.0000    0.2924   3.1416
  0.5    upwind      13.17    0.0000    0.1954   3.1416
  0.5      SUPG       0.00   -0.0000    0.2704   3.1416
  1.0  Galerkin      13.53    0.0000    0.5774   3.1416
  1.0    upwind      19.80    0.0000    0.2924   3.1416
  1.0      SUPG       0.00   -0.0000    0.4428   3.1416
  2.0  Galerkin      35.17   -0.3334    1.1010   3.1416
  2.0    upwind      18.17   -0.0000    0.3885   3.1416
  2.0      SUPG       0.00    0.0000    0.5572   3.1416
  5.0  Galerkin      69.61   -0.6961    2.1517   3.1416
  5.0    upwind       9.09    0.0000    0.4836   3.1416
  5.0      SUPG       0.00    0.0000    0.5773   3.1416
 
vel = 0 : skew/sym = 0.00e+00,  loop W = -4.29e-16
pi * u  = 3.1416

Pe_h 슬라이더를 1 너머로 밀면 오른쪽 a_E가 붉게 넘어가며 절점이 0 아래로 가라앉는다.

Pe_h2.0beta0.00

Leave beta at 0 and drag Pe_h past 1: a_E flips sign, and the marching profile starts ringing — at Pe_h = 2 the node next to the outlet dives to 0.000 with the boundary values still pinned at 0 and 1. Hit SUPG and the error goes to 0.00 %. The pink line at bottom right is the skew part of the matrix: none of the three buttons moves it.

beta를 0에 둔 채 Pe_h가 2가 되면 표의 min phi처럼 -0.3334 근처가 보인다. SUPG를 누르면 오차는 0 %로 붙는데, 분홍 반대칭 줄은 3.1416에서 움직이지 않는다.

경계값은 0과 1로 고정되어 있는데 $Pe_h = 2$ 의 갈러킨 해는 $-0.3334$ 까지 내려간다. $Pe_h = 5$ 에서는 $-0.6961$ 이다. 물리적으로 나올 수 없는 값이다.

$Pe_h = 1$ 에서 계수 하나가 부호를 바꾼다#

격자 Péclet 수를 $Pe_h = u h / (2\epsilon)$ 로 정의하면 위 행의 오른쪽 계수는 이렇게 다시 쓰인다.

K_{i,i+1} = \frac{\epsilon}{h}\left(Pe_h - 1\right)

$Pe_h > 1$ 이면 비대각 성분이 양수가 된다. 이 순간 행렬은 M-행렬이기를 그만두고, 이산 최대원리 — 내부 해가 경계값이 만드는 구간을 벗어나지 않는다는 보증 — 도 함께 사라진다. 표에서 $Pe_h = 1.0$ 까지 min phi가 0이었다가 2.0에서 음수가 되는 이유가 이것이다. 위에서 a_E가 붉게 넘어간 순간이 바로 이 부호 변화다.

안정화는 대칭을 되돌리지 않는다#

여기서 흔한 오독이 있다. 상류차분이나 SUPG가 "안정성을 회복시킨다"는 말은, 대칭성과 최소화 원리를 되돌린다는 뜻이 아니다. 표의 loop W 열을 보자. 세 스킴, 네 가지 $Pe_h$ — 총 12행 전부 3.1416이다. 이 값은 $\pi u$ 이고, $\epsilon$ 도 $\beta$ 도 들어 있지 않다.

이유는 간단하다. $\beta$ 가 키우는 것은 $\epsilon_{\text{eff}}$ 뿐이고, 그건 행렬의 대칭 부분이다. 반대칭 부분 $u$ 는 손대지 않는다. skew/sym 비가 $Pe_h = 5$ 에서 2.1517에서 0.4836으로 줄어든 것도 반대칭이 작아져서가 아니라 분모인 대칭 부분이 커져서다.

$Pe_h = 0.5$ 에서 갈러킨 max|err|%는 3.45, min phi는 0.0000, skew/sym은 0.2924다. $Pe_h = 1.0$ 이면 오차는 13.53까지 오르지만 최소값은 아직 0.0000이다. 식 $K_{i,i+1} = (\epsilon/h)(Pe_h - 1)$ 에서 $Pe_h = 1$ 이면 오른쪽 계수가 0이고, 1을 넘는 순간 양수가 된다. $Pe_h = 2$ 에서 갈러킨의 최소값은 -0.3334, 오차는 35.17 %, skew/sym은 1.1010이다. 같은 줄의 상류차분은 최소값 -0.0000, 오차 18.17 %, skew/sym 0.3885다. SUPG는 오차 0.00, skew/sym 0.5572. 세 줄의 loop W는 전부 3.1416이다. $u = 1$ 일 때 $\pi u$ 가 3.1416이라, $\beta$ 가 0이든 1이든 $\coth(Pe_h) - 1/Pe_h$ 이든 고리 일은 같다. $Pe_h = 5$ 의 갈러킨은 최소값 -0.6961, 오차 69.61 %, skew/sym 2.1517이다. 상류차분은 오차를 9.09 %로 내리고 skew/sym을 0.4836으로 줄인다. 줄어든 것은 분모의 대칭 부분이다. 속도 0인 대조는 skew/sym 0.00e+00, loop W -4.29e-16이다. 요소 10개, $u = 1$ 이면 $h = 0.1$ 이고 $Pe_h = uh/(2\epsilon)$ 라서 $Pe_h = 2$ 의 $\epsilon$ 은 0.025다. 그 $\epsilon$ 을 $\beta$ 로 키워도 $K_{i,i+1} - K_{i+1,i} = u$ 는 1로 남는다.

그래서 안정화가 실제로 사는 것은 더 약한 보증이다. 최소화 원리(에너지 노름에서의 최적 근사)를 잃은 대신, M-행렬 성질과 이산 최대원리를 산다. 값은 정확도로 치른다. $Pe_h = 2$ 에서 상류차분은 진동을 없앴지만 오차는 18.17 %다. 반면 SUPG는 $\beta$ 를 딱 필요한 만큼만 넣어서 절점에서 정확하다.

\beta_{\text{opt}} = \coth(Pe_h) - \frac{1}{Pe_h}

이 값은 $Pe_h \to 0$ 에서 0으로, $Pe_h \to \infty$ 에서 1로 간다. 확산이 지배하면 안정화를 끄고, 이류가 지배하면 완전 상류차분으로 간다는 뜻이다. 다만 절점 정확성은 1차원 상수계수 문제의 특권이다. 2차원에서는 유선 방향으로만 인공확산을 넣는 SUPG의 원래 형태가 필요하고, 그때도 "정확"은 사라진다.

유한체적법이 같은 자리에서 쓰는 다른 이름#

이 계산은 유한요소법 언어로 했지만 결론은 이산화 방법과 무관하다. 유한체적법에서 중심차분 대류항은 정확히 같은 스텐실을 만들고, 같은 $Pe_h = 1$ 에서 계수가 부호를 바꾼다. 그래서 등장하는 것이 1차 상류차분이고, 그것이 넣는 수치확산은 $\beta = 1$ 일 때의 $u h / 2$ 와 같은 값이다.

두 세계가 갈라지는 지점은 보존성이다. 유한체적법의 상류차분은 면 플럭스를 바꾸는 방식이라 총량 보존이 유지된다. 보존형과 원시형이 충격파 속도에서 갈라지던 것과 같은 종류의 분기다. 유한요소법의 인공확산은 강성행렬에 항을 더하는 것이라, 무엇을 보존하는지 따로 따져야 한다.

가중잔차 계열 안에서 다른 선택지를 고르면 또 다른 대가를 치른다. 최소자승법은 잔차 노름을 최소화하므로 항상 대칭 양정치 행렬을 만든다 — 최소화 원리가 돌아온다. 대신 조건수가 제곱되고, 선형 요소에서는 요소 내부의 2차 미분이 0이라 확산항이 통째로 사라진다. 기저와 적분 규칙이 정식화의 성질을 조용히 바꾸는 자리다.

$Pe_h$ 가 2면 조밀화 전에 유계성을 산다#

강의노트의 "어떤 방법을 써도 같다"는 문장은 자기수반 연산자 안에서만 맞다. 확산, 탄성, 포텐셜 유동은 그 안이다. 이류가 들어오면 $K_{i,i+1} - K_{i+1,i} = u$ 라서 대칭이 깨지고, loop W는 $\epsilon$ 을 어떻게 바꿔도 3.1416, 곧 $\pi u$ 에 고정된다. 그때 CG가 약속하던 에너지 노름의 최적 근사는 없다. 비대칭 행렬은 Krylov 쪽, Jacobi에서 BiCGSTAB까지로 넘기는 편이 정직하다.

min phi가 음수면 셀 Péclet가 1을 넘은 것이다. $Pe_h = 2$ 에서 갈러킨은 -0.3334, 오차 35.17 %다. 상류차분은 진동을 지우지만 오차 18.17 %를 남긴다. SUPG는 절점 오차 0.00인데 skew/sym은 1.1010에서 0.5572로 줄었을 뿐 반대칭 자체는 그대로다. 3차원에서 $Pe_h$ 를 5에서 1로 내리려면 셀이 125배다. 그 배수를 감당할 수 없을 때 사는 것은 정확도가 아니라 유계성이다.