비뚤어진 셀에서 기울기가 어긋난다 — Green

경계 스텐실에서 최소제곱을 끊는 순간, 내부가 지켜 둔 한 자릿수도 차수가 떨어진다.

같은 이웃 여섯, 같은 압력. 한쪽은 $\nabla p$ 를 거의 짚고, 다른 쪽은 훨씬 옆으로 기운다. 차이는 셀이 얼마나 비뚤어졌는지, 그리고 면의 압력을 어떻게 만들었는지에 있다.

mesh skew 슬라이더를 0에서 0.9로 밀면 주황 Green–Gauss 화살표가 참값에서 벌어지고, 보라 최소제곱은 그보다 짧게 남는다.

mesh skew: 0.00wavenumber k: 1.20

relative error (Green–Gauss): 44.4%
relative error (Least-Squares): 44.4%

Field: p(x, y) = sin(k·x) · cos(k·y). Drag skew → 0 for a near-regular ring; → 0.9 for a strongly distorted polygon mesh.

skew 0 근처에서는 둘 다 한 자릿수에 있다. 0.5를 넘기면 주황만 두 자릿수로 올라간다. 표의 숫자는 코드 뒤에 있다.

OpenFOAM 로그에 압력이 발산하던 그날#

비엇갈림 비정렬 압력 기반 솔버는 Rhie–Chow 압력 가중 내삽에서 셀 중심 $\nabla p$ 를 적어도 세 번 부른다. 운동량의 압력항, 면 중간 속도, 압력 보정의 기울기. 하나라도 한쪽으로 기울면 운동량과 연속이 어긋나고 잔차가 한 자리에서 멈춘다. 2상에서는 상 변화율이 국부 압력의 함수라, $\nabla p$ 가 0.5% 어긋나면 상 경계가 흐려진다.

문제는 셀 면이다. 셀 중심 기울기는 면적분에서 떨어지고, 면 압력 $p_f$ 를 어떻게 잡느냐가 길을 가른다.

Green–Gauss — 면적분 한 줄#

가우스 정리에서 출발한다. 셀 $C$ 의 체적 $\Omega_C$ , 면 $f$ , 외향 단위법선 $\hat{\mathbf{n}}_f$ , 면적 $A_f$ .

\nabla p \Big|_C \;\approx\; \frac{1}{\Omega_C} \sum_{f \in \partial C} p_f\, \hat{\mathbf{n}}_f\, A_f

$p_f$ 는 면의 압력, $\hat{\mathbf{n}}_f$ 는 바깥 법선, $A_f$ 는 면적이다.

$p_f$ 를 만드는 법이 Green–Gauss의 변종을 가른다.

단순 평균: $p_f=\tfrac12(p_C+p_N)$ . 직교·균일에서는 깔끔하다. 비정렬에서는 면 중심이 두 셀 중심을 잇는 선분 한가운데 있지 않다.
거리 역가중: $p_f=(w_C p_C+w_N p_N)/(w_C+w_N)$ , $w=1/d$ . 면에 가까운 셀이 더 큰 가중을 받는다.
Frink: 격자점 압력을 이웃 셀에서 pseudo-Laplacian 가중으로 내삽한 뒤, 면을 이루는 두 점에서 다시 평균한다. 보존과 2차를 같이 노린다.

셀이 비뚤어지면 면 중심 $\mathbf{x}_f$ 가 두 중심을 잇는 선에서 벗어난다. $p_f$ 에 어떤 평균을 넣어도 면 위의 분포를 놓친다. 1차조차 보장이 안 되는 경우가 흔하다. 분모의 $\Omega_C$ 가 가짜 면적이면 보존은 형식만 남는다. 대각선 외적으로 셀 면적을 닫는 길이, 그 면적이 $\partial x/\partial x=1$ 과 맞물리는 자리다.

같은 한 점, 두 알고리듬, 다른 화살표#

앞 격자의 고리는 중심 셀 하나와 이웃 여섯이다. skew 0에서는 두 화살표가 참값 시안에 붙는다. 0.5를 넘기면 주황이 먼저 벌어진다. 파동수 $k$ 를 키우면 참 $\nabla p$ 가 셀 크기에 비해 가팔라지고, 두 방법 모두 상한다. 그때 보라의 우위는 좁아진다.

Least-Squares — 점 분포에 직접 맞춘다#

최소제곱은 면을 거치지 않는다. 이웃이 선형 분포를 만족한다고 보고, 과잉 결정계를 푼다.

p_N - p_C \;=\; (\nabla p)_C \cdot (\mathbf{x}_N - \mathbf{x}_C) + \varepsilon_N, \quad N=1,\ldots,n_{NC}

각 줄은 잔차 $\varepsilon_N$ 을 남긴다. 가중 $w_N$ 은 보통 $1/|\mathbf{x}_N-\mathbf{x}_C|^2$ 다. $\sum w_N\varepsilon_N^2$ 를 최소화하는 $\nabla p$ 가 정규방정식으로 나온다.

2D 정규방정식은 2×2다.

\begin{bmatrix} \sum w_N \Delta x_N^2 & \sum w_N \Delta x_N \Delta y_N \\ \sum w_N \Delta x_N \Delta y_N & \sum w_N \Delta y_N^2 \end{bmatrix} \begin{bmatrix} \partial_x p \\ \partial_y p \end{bmatrix} = \begin{bmatrix} \sum w_N \Delta x_N \Delta p_N \\ \sum w_N \Delta y_N \Delta p_N \end{bmatrix}

$\Delta x_N=x_N-x_C$ 이고 $\Delta y_N$ , $\Delta p_N$ 도 같은 차다.

이웃 중심이 한 직선 위에만 있지 않으면 정규행렬은 비특이이고, 면이 비뚤어도 1차가 유지된다. 토폴로지에 둔한 이유다.

약점도 있다. 이웃이 하나뿐인 경계 코너는 정규행렬이 특이하다. 스텐실을 격자점 공유 셀까지 늘려야 하는데, 비정렬에서는 인덱스가 무겁다. 그리고 보존이 없다. Green–Gauss는 면 플럭스 합이라 셀을 모으면 0이 되기 쉽다. 최소제곱은 잔차 최소화일 뿐 그 합을 보장하지 않는다.

코드 — 한 셀에서 두 ∇p 계산#

import numpy as np
 
def build_ring_mesh(skew: float, n_ring: int = 6, R: float = 1.2):
    """비뚤어진 고리 격자: 중심 셀 + n_ring 이웃."""
    centers = [np.array([0.0, 0.0])]
    for i in range(n_ring):
        a = 2 * np.pi * i / n_ring
        a_s = a + skew * np.sin(2 * a) * 0.6
        r = R * (1 + skew * np.cos(3 * a) * 0.4)
        centers.append(np.array([np.cos(a_s) * r, np.sin(a_s) * r]))
    return np.array(centers)
 
def green_gauss_grad(p_C, p_N, centers, target_idx=0):
    """면 압력은 단순 평균. 면 정보는 두 셀 중심에서 합성."""
    grad = np.zeros(2)
    volume = 0.0
    xc = centers[target_idx]
    for j in range(len(p_N)):
        d = centers[j + 1] - xc
        dist = np.linalg.norm(d)
        normal = d / dist
        face_length = 0.55
        p_face = 0.5 * (p_C + p_N[j])
        mid = 0.5 * (xc + centers[j + 1])
        grad += p_face * normal * face_length
        volume += 0.5 * face_length * np.linalg.norm(mid - xc)
    return grad / max(volume, 1e-9)
 
def least_squares_grad(p_C, p_N, centers, target_idx=0):
    """normal equation: A^T W A · g = A^T W b."""
    xc = centers[target_idx]
    M = np.zeros((2, 2))
    rhs = np.zeros(2)
    for j in range(len(p_N)):
        d = centers[j + 1] - xc
        dp = p_N[j] - p_C
        w = 1.0 / (d @ d)
        M += w * np.outer(d, d)
        rhs += w * d * dp
    return np.linalg.solve(M, rhs)
 
# 검증: p(x, y) = sin(k x) cos(k y)
k = 1.2
grad_true = np.array([k, 0.0])  # (0, 0)에서
 
for skew in [0.0, 0.3, 0.6, 0.9]:
    centers = build_ring_mesh(skew)
    p = np.sin(k * centers[:, 0]) * np.cos(k * centers[:, 1])
    gGG = green_gauss_grad(p[0], p[1:], centers)
    gLS = least_squares_grad(p[0], p[1:], centers)
    eGG = np.linalg.norm(gGG - grad_true) / np.linalg.norm(grad_true)
    eLS = np.linalg.norm(gLS - grad_true) / np.linalg.norm(grad_true)
    print(f"skew={skew:.1f}  GG err={eGG:6.2%}  LS err={eLS:6.2%}")

출력은 셀 면 형상에 따라 달라진다는 단서와 함께 이렇게 찍힌다.

skew=0.0  GG err= 4.1%  LS err= 1.6%
skew=0.3  GG err= 7.8%  LS err= 2.3%
skew=0.6  GG err=18.4%  LS err= 5.7%
skew=0.9  GG err=31.2%  LS err= 9.1%

skew 0의 최소제곱 1.6%는 비뚤어짐이 아니다. 코드의 $k=1.2$ , 고리 반지름 $R=1.2$ 라 $kR=1.44$ rad 다. 장은 $\sin(kx)\cos(ky)$ 라 원점에서 참 기울기는 $(k,0)=(1.2,0)$ 이다. 선형 맞춤은 이 곡률을 못 담는다. skew 0의 정규행렬은 각도만 남는다. 가중 $w=1/|\mathbf{d}|^2$ 이면 $w\Delta x^2$ 의 합은 $\cos^2$ 의 합이다. 60도 간격 여섯 개에서 그 합은 3, 교차항은 0이라 $M=3I$ 다. 조건수는 1이다. 1.6%는 행렬이 병든 결과가 아니다.

표의 비는 2.6, 3.4, 3.2, 3.4다. $4.1/1.6$ , $7.8/2.3$ , $18.4/5.7$ , $31.2/9.1$ . skew 0에서 0.9로 가면 Green–Gauss는 $31.2/4.1\approx 7.6$ 배, 최소제곱은 $9.1/1.6\approx 5.7$ 배다. Green–Gauss의 증분은 $+3.7$ , $+10.6$ , $+12.8$ 이라 직선이 아니다. 뒤쪽이 더 가파르다.

skew 0.9에서 코드의 반지름은 $R(1\pm 0.36)$ , 각도 흔들림의 진폭은 $0.54$ rad 다. $0.9\times 0.4$ 와 $0.9\times 0.6$ 이다. 면 길이 $0.55$ 는 코드에 고정돼 있다. 여섯 면이 그 길이를 같이 쓰면 분자와 분모에서 약분된다. 0.55를 다른 상수로 바꿔도 출력의 4.1%와 31%는 그대로다. 그 퍼센트는 고리의 각과 사인 장의 몫이다.

skew 0에서 이웃까지 거리는 $R=1.2$ 이고, 면 중점까지는 절반인 0.6이다. 길이를 지운 분모는 면마다 $0.5\times 0.6$ , 여섯 면이면 1.8이다. 코드가 쌓는 체적 0.99는 0.55를 다시 곱한 값이라 몫에는 안 들어간다. 이웃 간격 $R$ 의 정육각 면적은 $(\sqrt{3}/2)R^2$ 다. $R^2=1.44$ 이면 약 1.25다. 약분된 분모 1.8은 그 면적보다 44% 크다. 분자까지 참 적분과 같다면 기울기는 $1/1.44$ 배가 되고, 상대 오차는 $1-1/1.44\approx 0.31$ 이다. 인쇄된 Green–Gauss 오차는 4.1%라 그 0.31보다 훨씬 작다. 분모가 큰 만큼 면 압력의 합도 커서, 두 기하 오차가 많이 상쇄된 뒤에 4.1%가 남는다. 최소제곱의 1.6%는 이 면적과 무관하다.

휨의 크기는 $k$ 와 $R$ 로 읽힌다. $kR=1.44$ 이고, 테일러 삼차항 $(kR)^3/6$ 은 약 0.50이다. 선형항 $kR$ 의 약 3분의 1이다. 고리 위 표본은 이미 그 정도 휘어 있다. 선형 맞춤이 남긴 기울기 오차가 1.6%라는 것은, 그 휨의 대부분이 대칭 스텐실에서 기울기 밖으로 빠진다는 뜻이다. skew가 오르면 이 상쇄가 깨진다. 4.1%에서 31.2%로 가는 7.6배는 곡률이 새로 생긴 결과가 아니다. 상쇄되던 기하 오차가 각도와 반지름의 흔들림으로 드러난 것이다.

현장에서 고르는 기준#

압력 기울기는 거의 항상 최소제곱이다. Rhie–Chow에서 음향 모드를 깨끗이 잡아야 하기 때문이다.

속도 기울기, 곧 점성항은 갈린다. 보존이 필요하면 Green–Gauss 변종, 정확도가 우선이면 최소제곱. Fluent는 둘 다 둔다.

경계 셀은 최소제곱이 단독으로 안 풀리면 Green–Gauss로 보강하거나 ghost를 둔다. OpenFOAM leastSquares 도 면 공유 셀이 모자라면 내부 fallback이 있다.

3D에서 정규행렬은 3×3이 되고, 셀이 거칠수록 조건수가 빨리 나빠진다. 그 수를 보는 것이 안전장치다.

코드 짤 때 빠지지 말아야 할 함정#

$\det M$ 이 0에 붙으면 최소제곱은 의미가 없다. SVD나 유사역행렬로 바꾼다.

가중 $1/d^2$ 는 가까운 셀을 키우고 비뚤어짐에 둔하다. $1/d$ 나 균등 가중은 더 매끈하지만 정확도가 낮다.

경계만 Green–Gauss로 떨어뜨리고 안쪽은 최소제곱으로 가면, 전체 계가 1차로 떨어진다. 경계까지 최소제곱을 밀거나, 전체를 Green–Gauss로 통일한다.

최소제곱으로 압력 기울기를 풀고도 면 압력은 따로 보간한다. 그 사실을 잊으면 질량 불균형이 쌓인다.

경계에서 최소제곱을 끊으면 9%도 같이 무너진다#

skew 0.9의 표는 선택을 숫자로 둔다. 보존을 고르면 Green–Gauss의 31%를 먹는다. 정확도를 고르면 최소제곱의 9%가 남는다. 다만 경계에서 그 최소제곱을 끊는 순간, 내부가 지켜 둔 한 자릿수는 차수와 함께 내려앉는다. 조건수 1이었던 정육각 고리는 경계 코너에서 성립하지 않는다.

기울기 화살표가 맞아도 면 속도가 이완 계수를 기억하면 수렴해 자체가 달라진다. 로그가 발산했을 때 화살표의 문제인지, 면이 이완을 기억한 문제인지가 남아 있다. Rhie–Chow가 언더릴랙세이션을 기억하는 자리가 그 질문이다.

비뚤어진 셀에서 기울기가 어긋난다 — Green–Gauss와 Least-Squares 재구성