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

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

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

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

그 문장에는 조건이 하나 붙어 있는데, 막대 문제 안에 숨어 있어서 보이지 않는다. 유동 방정식에는 그 조건이 없다. 이 글에서는 그 조건이 무엇이고, 조건이 깨졌을 때 정확히 어떤 보증을 잃는지를 잰다. 그리고 안정화 기법이 되돌려 주는 것이 무엇인지 — 그리고 무엇은 절대 되돌려 주지 못하는지 — 를 숫자로 확인한다.

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

최소 포텐셜 에너지 정식화는 총 포텐셜 에너지를 절점 변위로 미분해서 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이다.

아래 시뮬레이션에서 자유도 2개짜리 장난감 행렬로 직접 확인해보자.

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.

advection a를 0에 두면 분홍 점이 한 바퀴 돌고 와서 0.000을 들고 온다. 회색 타원은 등에너지선이고, 노란 공은 그 선을 가로지르지 않고 곧장 안쪽으로 미끄러진다. a를 올리면 한 바퀴의 일이 정확히 $2\pi a$ 가 되고, 공은 나선을 그리며 등에너지선을 거슬러 올라가기 시작한다. 이 상태에서 diffusion s를 끝까지 밀어보자. 한 바퀴의 일은 꿈쩍도 하지 않는다.

이류항이 깨뜨리는 양은 정확히 $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 행렬로 접을 수 있었던 것과는 사정이 다르다.

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

경계값은 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에서 음수가 되는 이유가 이것이다.

아래 시뮬레이션은 같은 계를 시간 전진으로 푼다. 정상해가 완성되는 과정을 보면 진동이 어디서 자라나는지가 보인다.

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를 1 너머로 밀어보자. 오른쪽 a_E 막대가 왼쪽으로 넘어가면서 빨갛게 바뀌고, 그 다음 시간 전진에서 절점 값이 0 아래로 내려간다. SUPG 버튼을 누르면 오차가 0 %가 된다. 그런데 오른쪽 아래 분홍 줄 — 행렬의 반대칭 성분 — 은 세 버튼 어느 것으로도 움직이지 않는다.

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

여기서 흔한 오해를 정리해야 한다. 상류차분이나 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으로 줄어든 것도 반대칭이 작아져서가 아니라 분모인 대칭 부분이 커져서다.

그래서 안정화가 실제로 사는 것은 더 약한 보증이다. 최소화 원리(에너지 노름에서의 최적 근사)를 잃은 대신, 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이라 확산항이 통째로 사라진다. 불연속 갈러킨의 적분점 하한을 따지던 글에서처럼, 기저와 적분 규칙이 정식화의 성질을 조용히 바꾸는 자리다.

비대칭 행렬을 넘겨받았을 때#

강의노트의 "어떤 방법을 써도 같다"는 문장의 유효 범위는 자기수반 연산자다. 확산, 탄성, 포텐셜 유동은 그 안에 있다. 이류가 들어오는 순간 밖으로 나간다.

실무에서는 세 가지를 순서대로 확인하면 된다. 첫째, 조립된 행렬이 대칭인가. 대칭이면 CG 계열을 쓸 수 있고 에너지 노름에서 최적 근사라는 보증이 따라온다. 비대칭이면 GMRES 계열이고, 그 보증은 없다. 둘째, 셀 Péclet 수가 1을 넘는가. 넘으면 진동은 버그가 아니라 스킴의 정의된 거동이다. 셋째, 안정화를 켰다면 그것이 산 것이 정확도인지 유계성인지. 대개는 유계성이고, 정확도는 지불한 쪽이다.

해가 경계값 구간 밖으로 나갔을 때 격자를 조밀하게 하는 것은 임시방편이 아니라 정공법이다. $h$ 를 줄이면 $Pe_h$ 가 같이 줄기 때문이다. 다만 3차원에서 $Pe_h$ 를 5에서 1로 내리려면 셀 수가 125배가 된다. 그 계산을 해본 다음에 안정화 항을 고르는 편이 낫다.