격자가 비뚤어지면 라플라시안이 새어 나간다

직교 보정의 수렴 벽은 70°가 아니라, cosθ가 1/2이 되는 60°에 있다.

정육면체에서 잘 돌던 확산 솔버를 비뚤어진 삼각형으로 옮기면, 답이 슬그머니 틀린다. 방정식은 같고 코드도 같다. 온도장만 흐려지고 수렴은 느려진다. 경계조건도 아니고 물성도 아니다. 면 위에서 기울기를 재는 방식이다.

비직교각 슬라이더를 20°에서 65°로 밀면 과완화의 청록 벡터는 면적보다 길어지고, 최소 보정의 분홍 벡터는 더 길어진다.

non-orthogonality θ = 35°

over-relaxed   |E_f|/|S_f| = 1.221   |T_f|/|S_f| = 0.700
explicit/implicit ratio ρ = |T_f|/|E_f| :   min 0.70  |  orth 0.60  |  over 0.57

청록이 분홍보다 짧아지면 화면 아래 ρ가 1을 넘는다. 45°에서 최소 보정만 먼저 빨간 쪽으로 가는지 보면 된다.

면 위의 기울기가 말썽이다#

유한체적법(FVM, 셀 단위로 보존량을 적분하는 방법)에서 확산항은 면을 지나는 플럭스의 합이다. 면 $f$ 하나의 확산 플럭스는 이렇게 쓴다.

F_f = \Gamma_f \, \nabla\phi_f \cdot \mathbf{S}_f

$\Gamma_f$ 는 면의 확산계수다. $\nabla\phi_f$ 는 면 기울기다. $\mathbf{S}_f$ 는 바깥을 향한 면적벡터, 넓이 곱하기 법선이다.

$\nabla\phi_f \cdot \mathbf{S}_f$ 를 셀 값으로 내리는 일이 문제다. 정렬 직교격자에서는 $\mathbf{d}$ , 두 셀 중심을 잇는 선이 법선과 나란하다. $(\phi_N-\phi_P)/|\mathbf{d}|$ 가 곧 법선 기울기다.

비정렬에서는 그 가정이 깨진다. $\mathbf{d}$ 와 $\mathbf{S}_f$ 가 어긋난다. 그 각이 비직교각 $\theta$ 다. $(\phi_N-\phi_P)/|\mathbf{d}|$ 는 이제 $\mathbf{d}$ 방향 기울기지, 법선 기울기가 아니다. 차이를 무시하면 확산이 새고, 해는 흐려진다.

S_f를 두 조각으로 쪼갠다#

면적벡터를 두 성분으로 나눈다. 하나는 $\mathbf{d}$ 방향, 하나는 나머지가다.

\mathbf{S}_f = \mathbf{E}_f + \mathbf{T}_f

$\mathbf{E}_f$ 는 $\mathbf{d}$ 와 나란한 성분이다. $\mathbf{T}_f$ 는 접선 쪽 나머지다. 플럭스에 넣으면 항이 둘로 갈린다.

F_f = \Gamma_f |\mathbf{E}_f| \frac{\phi_N - \phi_P}{|\mathbf{d}|} + \Gamma_f \, \overline{\nabla\phi}_f \cdot \mathbf{T}_f

앞항은 $\phi_N-\phi_P$ 만 쓴다. 음함수로 행렬에 넣는다. 셀 중심 두 점뿐이다. 뒷항은 면 평균 기울기 $\overline{\nabla\phi}_f$ 가 필요하다. 이웃 기울기를 평균해 얻는다. 계수로 넣기 어려워 명시적 보정으로 우변에 던진다. 그 기울기를 최소제곱으로 다시 짜는 쪽은 비정렬에서 기울기를 살리는 이동 최소제곱에 따로 있다. 지연 보정이 흔들릴 때 범인이 ρ가 아니라 그 기울기인 경우가 있다.

$\theta=0$ 이면 $\mathbf{T}_f=0$ 이다. 옛 중심차분으로 돌아간다.

최소·직교·과완화: 세 가지 분해#

$\mathbf{E}_f$ 를 $\mathbf{d}$ 방향으로 두는 길이는 하나가 아니다. $\hat{\mathbf{d}}=\mathbf{d}/|\mathbf{d}|$ , $\cos\theta=\hat{\mathbf{d}}\cdot\mathbf{S}_f/|\mathbf{S}_f|$ 로 두면 세 가지다.

최소 보정: $|\mathbf{E}_f|=|\mathbf{S}_f|\cos\theta$ . $\mathbf{T}_f$ 가 가장 짧다.
직교 보정: $|\mathbf{E}_f|=|\mathbf{S}_f|$ . 음함수 크기를 면적 그대로 둔다.
과완화: $|\mathbf{E}_f|=|\mathbf{S}_f|/\cos\theta$ . 비직교가 커질수록 음함수 성분을 더 키운다.

앞 화면에서 θ를 60°까지 올리면 과완화의 $\mathbf{E}_f$ 가 면적벡터보다 길어진다. 대각 계수가 커진다는 신호다. 최소 보정은 반대로 $\mathbf{E}_f$ 가 줄고 $\mathbf{T}_f$ 가 길어진다.

코드: 지연 보정 라플라스 솔버#

명시적 보정은 직전 반복의 기울기를 쓴다. 한 번 풀고 끝이 아니다. 보정을 갱신하며 여러 번 쓸어 주는 지연 보정 루프다.

import numpy as np
 
def decompose_face(S_f, d, scheme="over"):
    """면적벡터 S_f = E_f(음함수, d 방향) + T_f(명시적 보정)."""
    d_hat = d / np.linalg.norm(d)
    Smag = np.linalg.norm(S_f)
    dS = d_hat @ S_f                      # = |S_f| cos(비직교각)
    if scheme == "min":                   # 최소 보정
        E = dS * d_hat
    elif scheme == "orth":                # 직교 보정
        E = Smag * d_hat
    else:                                 # 과완화 (over-relaxed)
        E = (Smag**2 / dS) * d_hat
    return E, S_f - E                     # E_f, T_f
 
def sweep_residual(S_f, d, scheme, sweeps=25):
    """E_f는 음함수, T_f는 지연. 잔차는 ρ=|T|/|E| 비율로 수축한다."""
    E, T = decompose_face(S_f, d, scheme)
    rho = np.linalg.norm(T) / np.linalg.norm(E)
    r, hist = 1.0, [1.0]
    for _ in range(sweeps):
        r *= rho                          # 지연 보정 1회 = 잔차 × ρ
        hist.append(r)
    return rho, hist
 
S_f = np.array([1.0, 0.0])                # 면 법선 (면적벡터, |S|=1)
for ang in (20, 45, 65):
    d = np.array([np.cos(np.radians(ang)), np.sin(np.radians(ang))])
    for s in ("min", "orth", "over"):
        rho, hist = sweep_residual(S_f, d, s)
        tag = "수렴" if rho < 1 else "발산"
        print(f"θ={ang:2d}°  {s:5s}  ρ={rho:.3f}  {tag}  (25회 후 {hist[-1]:.1e})")

출력은 각이 커질수록 세 방식이 갈라진다는 뜻이다. 20°에서는 셋 다 무난하다. 65°에서 최소 보정은 $\rho=\tan 65°\approx 2.14$ 로 발산하고, 과완화는 $\rho=\sin 65°\approx 0.91$ 로 아직 수렴한다. 같은 격자, 같은 방정식인데 분해 하나가 수렴과 발산을 가른다.

세 분해의 ρ는 각도로 닫힌다. 면적 길이를 1로 두면 최소 보정은 E의 길이가 cosθ, T의 길이가 sinθ라 ρ는 tanθ다. tan 45°는 1이다. 45°가 수렴의 벽이다. 그 각에서 지연 루프는 잔차에 1을 곱한다. 25번을 돌려도 숫자가 그대로다. 20°의 tan은 약 0.364다. 65°의 tan은 약 2.14다. 발산 쪽 숫자는 이 탄젠트에서 나온다.

직교 보정은 E의 길이가 1이고 T의 길이는 2에서 2cosθ를 뺀 값의 제곱근이다. ρ가 1이 되려면 cosθ가 1/2, 즉 60°다. 65°에서는 이미 약 1.07이다. 70°까지 가야 직교가 위험해진다는 말은 늦다. 벽은 60°에 있다.

과완화의 ρ는 sinθ다. 65°에서 약 0.906이다. 스니펫은 잔차에 이 비율을 25번 곱한다. 0.906의 25제곱은 약 0.085다. 25회 뒤에도 명시 오차의 약 8퍼센트가 남는다. 같은 식을 20°에 쓰면 sin이 약 0.342이고, 25제곱은 $10^{-12}$ 수준으로 꺼진다. 보정 횟수를 2나 3으로 두는 습관은 65° 격자에서 0.906의 세제곱, 약 0.74다. 오차의 약 4분의 3이 그 시간 스텝에 남는다.

코드의 각은 d와 S 사이의 각이다. S를 단위 벡터로 고정하면 20°, 45°, 65°가 비교점이다. 최소 보정의 45°는 tan이 정확히 1이라, 경험으로 고른 문턱이 아니다.

왜 ρ가 T 길이 나누기 E 길이인가. 지연 보정은 음함수 계수, E 길이로 풀고, 매 쓸기마다 T 길이만큼의 오차를 명시적으로 되먹인다. 한 번 쓸 때마다 잔차가 그 비율로 줄거나 는다. 비율이 1보다 작아야 루프가 닫힌다.

비직교가 커지면 무슨 일이 일어나는가#

세 방식의 잔차를 한 화면에서 같이 민다. 각 슬라이더가 ρ를 바꾸고, 쓸기마다 곡선이 그 비율로 내려앉는다.

mesh non-orthogonality θ = 55°

● minimum ρ=1.43 diverges● orthogonal ρ=0.92 converges● over-relaxed ρ=0.82 converges

30° 아래에서는 셋이 같이 떨어진다. 45°에서 최소 보정만 수평이 된다. 60°를 넘기면 직교도 ρ가 1을 넘고, 과완화만 sinθ라서 남는다.

대가는 있다. 과완화는 음함수 성분을 키우는 만큼 명시 보정도 크다. 같은 잔차에 닿으려면 쓸기를 여러 번 돌려야 한다. nNonOrthogonalCorrectors가 그 횟수다. 격자가 좋으면 0이나 1이면 된다. 심하게 찌그러지면 2나 3까지 올린다. 65°에서 그 두세 번은 오차의 4분의 3을 남겨 둔다. 무작정 키우면 한 시간 스텝이 비싸진다. 각 자체를 줄이는 편이 거의 항상 싸다.

45°에서 최소 보정을 버리는 이유#

45°에서 최소 보정은 잔차를 줄이지 못한다. ρ가 정확히 1이기 때문이다. 직교는 60°에서 같은 벽에 닿고, 65°에서는 이미 1을 넘는다. 과완화는 sinθ라 90° 직전까지 ρ가 1을 넘지 않는다. 대신 한 스텝에 남는 오차가 크다.

횟수를 더 줄 것인가, 각이 60°를 넘기 전에 격자를 다시 깔 것인가. 삼각형을 덮는 순서가 그 각을 만든다. 점 구름을 삼각형으로 덮는 두 생성법은 그 순서가 어디서 비직교를 남기는지 묻는다. 보정 루프를 늘리는 편이 싼가, 프론트를 다시 밀고 나가는 편이 싼가.

격자가 비뚤어지면 라플라시안이 새어 나간다 — 비직교 확산 플럭스 보정