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

정육면체 격자에서 잘 돌던 확산 솔버를 비뚤어진 삼각형 격자로 옮기면, 답이 슬그머니 틀어진다. 같은 방정식, 같은 코드인데 온도장이 흐릿해지고 수렴은 느려진다. 범인은 경계조건도 물성도 아니다. 면 위에서 기울기를 재는 방식이다. 이 글은 비정렬 격자에서 확산 플럭스를 어떻게 이산화하는지, 그리고 면적벡터를 두 조각으로 쪼개는 비직교 보정(non-orthogonal correction)이 왜 필요한지를 작은 솔버로 직접 만져보며 정리한다. 끝까지 읽으면 OpenFOAM 설정 파일에 늘 보이던 nNonOrthogonalCorrectors가 무엇을 세는 숫자인지 알게 된다.

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

유한체적법(FVM, 셀 단위로 보존량을 적분하는 방법)에서 확산항은 셀 표면을 지나는 플럭스의 합이다. 면 $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{E}_f$는 $\mathbf{d}$와 나란한 성분, $\mathbf{T}_f$는 그 나머지(접선 성분)다. 이 분해를 플럭스에 넣으면 항이 둘로 갈린다.

앞항은 $\phi_N - \phi_P$만 쓰므로 **음함수(implicit)**로 행렬에 넣는다. 셀 중심 두 점만 관여해 깔끔하다. 뒷항은 면 평균 기울기 $\overline{\nabla\phi}_f$가 필요하다. 이건 이웃 셀들의 기울기를 평균해 얻는다. 행렬 계수로 넣기 어려워 명시적(explicit) 보정으로 우변에 던진다.

핵심은 이것이다. 직교 성분 $\mathbf{E}_f$는 견고하게 음함수로 풀고, 비직교가 만든 군더더기 $\mathbf{T}_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|$로 두면:

• 최소 보정(minimum): $\mathbf{E}_f = (\hat{\mathbf{d}}\cdot\mathbf{S}_f)\,\hat{\mathbf{d}}$, 즉 $|\mathbf{E}_f| = |\mathbf{S}_f|\cos\theta$. $\mathbf{T}_f$ 길이가 가장 짧다.
• 직교 보정(orthogonal): $\mathbf{E}_f = |\mathbf{S}_f|\,\hat{\mathbf{d}}$. 음함수 성분 크기를 면적 그대로 유지한다.
• 과완화(over-relaxed): $\mathbf{E}_f = \dfrac{|\mathbf{S}_f|^2}{\hat{\mathbf{d}}\cdot\mathbf{S}_f}\,\hat{\mathbf{d}}$, 즉 $|\mathbf{E}_f| = |\mathbf{S}_f|/\cos\theta$. 비직교가 커질수록 음함수 성분을 더 키운다.

세 선택의 기하가 한눈에 들어오지 않는다. 아래 시뮬레이션에서 직접 조작해보자. 비직교각 슬라이더를 끌고, 세 방식을 번갈아 눌러 $\mathbf{E}_f$(청록, 음함수)와 $\mathbf{T}_f$(분홍, 보정)가 어떻게 달라지는지 보라.

$\theta$를 60°까지 올리면 과완화의 $\mathbf{E}_f$가 면적벡터보다 길어진다. 음함수 대각 계수가 커져 더 안정하다는 신호다. 반대로 최소 보정은 $\mathbf{E}_f$가 짧아지고 보정 $\mathbf{T}_f$가 급격히 길어진다. 화면 아래 ρ 값이 1을 넘으면(빨간색) 그 방식은 위험하다.

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

명시적 보정항은 직전 반복의 기울기를 쓴다. 그래서 한 번 풀고 끝이 아니라, 보정을 갱신하며 여러 번 쓸어주는 지연 보정(deferred correction) 루프가 된다. 면적벡터 분해와 그 수렴 거동을 작은 코드로 확인한다.

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$로 여전히 수렴한다. 같은 격자, 같은 방정식인데 분해 방식 하나가 수렴과 발산을 가른다.

왜 $\rho = |\mathbf{T}_f|/|\mathbf{E}_f|$일까. 지연 보정은 음함수 계수 $\propto |\mathbf{E}_f|$로 풀고, 매 쓸기마다 $\propto |\mathbf{T}_f|$만큼의 오차를 명시적으로 되먹인다. 한 번 쓸 때마다 잔차가 그 비율만큼 줄거나 는다. 이 비율이 1보다 작아야 루프가 닫힌다.

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

세 방식의 운명을 한 화면에 겹쳐 보자. 아래 시뮬레이션에서 직접 조작해보자. 격자 비직교각을 끌어 올리며 각 방식의 잔차가 몇 번의 쓸기 만에 떨어지는지 비교하라.

$\theta$가 30° 아래면 셋이 비슷하게 잘 떨어진다. 45°를 넘기면 최소 보정이 먼저 무너진다($\rho \ge 1$). 70°까지 밀어붙이면 직교 보정도 위태롭고, 과완화만 홀로 살아남는다. 과완화의 $\rho = \sin\theta$는 어떤 각에서도 1을 넘지 않기 때문이다. 이것이 상용·오픈소스 코드가 대부분 과완화를 기본값으로 쓰는 이유다.

대가는 있다. 과완화는 음함수 성분을 키우는 만큼 명시적 보정도 커서, 같은 잔차에 닿으려면 보정 쓸기를 여러 번 돌려야 한다. nNonOrthogonalCorrectors가 바로 이 쓸기 횟수다. 격자가 좋으면 01이면 충분하고, 심하게 찌그러진 격자에서는 23까지 올린다. 무작정 키우면 한 시간 스텝이 비싸지므로, 격자 품질을 먼저 손보는 편이 거의 항상 낫다.

기억할 점#

• 비정렬 격자의 확산항은 $\mathbf{S}_f = \mathbf{E}_f + \mathbf{T}_f$로 쪼개, $\mathbf{E}_f$는 음함수로 행렬에 넣고 $\mathbf{T}_f$는 지연 보정으로 우변에 던진다.
• 분해는 최소·직교·과완화 세 가지이며, 수렴은 $\rho = |\mathbf{T}_f|/|\mathbf{E}_f|$가 쥔다. 과완화의 $\rho=\sin\theta$는 늘 1 미만이라 가장 견고하다.
• 보정 쓸기 횟수(nNonOrthogonalCorrectors)를 늘리기 전에, 비직교각 자체를 줄이는 격자 개선이 우선이다.