[논문 리뷰] 폭풍이 오기 전에 끝내야 하는 시뮬레이션 — Casulli–Walters(2000) 비정렬 직교 격자 천수방정식

함부르크 BAW는 1990년 6월, Jade–Weser 강어귀 조위 13일을 다시 맞춰야 했다. 삼각형 15만 개. 시간 스텝은 5분. 실시간보다 빨리 끝나야 다음 예보에 손이 간다.

명시로 밀면 자유표면파 셀러리티 $c=\sqrt{gh}$ 가 CFL을 붙잡는다. 깊이 28 m에서 $c \approx 16.6$ m/s. 가장 짧은 변이 $\sqrt{10}$ m이면 허용 스텝은 0.2초다. 13일은 1,123,200초다. 0.2초로 자르면 5,616,000번, 논문의 560만 스텝과 같다. 5분은 3,744번. 비는 딱 1500배다. 300초를 0.2초로 나눈 값이라, 날짜를 바꿔도 이 배율은 그대로다.

그 배율을 버티는 쪽이 Casulli–Walters(2000)의 반음해다. 자유표면 구배만 새 시각에 두고, 이송은 옛 시각에 남긴다.

CFL 슬라이더를 1 위로 올리면 빨간 명시 곡선이 터지고, θ를 올리면 청록 곡선의 진폭만 깎인다.

CFL=6에서도 청록 봉우리는 남는다. θ를 1까지 밀면 그 봉우리가 빨리 죽고, 진폭이 버티는 실무 구간은 0.55에서 0.6이다.

논문 메타데이터#

• 저자: Vincenzo Casulli (Univ. of Trento), Roy A. Walters (USGS)
• 학술지: International Journal for Numerical Methods in Fluids, 32, 331–348 (2000)
• 제목: "An unstructured grid, three-dimensional model based on the shallow water equations"
• 키워드: semi-implicit, shallow water, unstructured orthogonal grid, wetting/drying

명시적 솔버가 강어귀를 못 푸는 이유#

깊이 10 m면 $c \approx 10$ m/s, 100 m면 $c \approx 31$ m/s. 유속 $u$ 는 보통 1 m/s 아래다. 조석이 강을 오가는 시간은 $L/u$ 인데, 명시 스텝은 $\Delta x/c$ 에 묶인다. $c/u \approx 10\sim30$. 물이 한 번 지나가는 동안 솔버는 시간을 그보다 잘게 썬다.

비정렬 격자는 그 제한을 가장 짧은 변에 몰아넣는다. Big Lost River는 1.2 km × 1 km에 삼각형 1.26만 개, 최단변 약 5 m. 명시 한계는 0.5초쯤이다. 5 m를 0.5초로 나누면 10 m/s. 깊은 골의 16 m/s가 시계를 멈추는 게 아니다. 얕은 곳의 짧은 변이 멈춘다.

느린 항은 명시, 빠른 항은 음해. 자유표면 압력 구배만 음해로 넘기면 $\sqrt{gh}$ 가 안정 조건에서 빠진다. 남는 제한은 유속 쪽이다.

직교 비정렬 격자 — Voronoi–Delaunay의 짝#

비정렬에 유한차분을 그대로 올리면 두 곳이 깨진다. 셀 간격이 모양마다 달라 구배가 한쪽으로 기운다. 곡선 좌표로 피하면 변환 항이 정확도와 안정성을 같이 깎는다.

Casulli와 Walters는 격자 쪽을 제한했다. 직교 비정렬 격자. 이웃 중심을 잇는 선분이 공유 변과 수직이다. Voronoi 다각형과 그 쌍대인 Delaunay 삼각분할이 이 조건을 만족한다. 예각 삼각형만일 때다. 둔각이 하나 들어가면 외심이 삼각형 밖으로 나가고, 중심 연결선은 변과 수직이어도 셀 안을 지나지 않는다. 빈 외접원을 지키는 점 삽입이 그 예각을 얼마나 오래 유지하는지가, 이 차분이 성립하는 시간이다.

직교가 주는 것은 차분 한 줄이다. 변의 수직 속도 $u_j$ 만 두면 $\partial \eta / \partial n$ 이 $(\eta_{i_2} - \eta_{i_1})/\delta_j$ 가 된다. 곡률 보정은 없다.

반음해 분리: θ 방법으로 두 시간 스케일 떼기#

θ 방법은 implicitness factor $\theta \in [1/2, 1]$ 로 시간을 가른다. 운동량에서는 자유표면 구배를 $\theta \eta^{n+1} + (1-\theta)\eta^n$ 으로, 자유표면식에서는 속도를 $\theta u^{n+1} + (1-\theta)u^n$ 으로 둔다.

이산화된 운동량이다. 연직 점성과 바람 응력이 들어 있다.

$\mathbf{U}_j^{n+1}$ 는 변 $j$ 의 수직 속도, $\mathbf{A}_j^n$ 은 연직 점성·바람·바닥 마찰을 모은 삼중대각, $\mathbf{G}_j^n$ 은 이송·코리올리·수평 마찰이다. 변마다 독립이라 $N_s$ 개의 $N_z \times N_z$ 삼중대각으로 갈라진다.

자유표면식은 이렇다.

$P_i$ 는 셀 면적, $\lambda_j$ 는 변 길이, $s_{i,l}$ 은 바깥 방향 부호다.

파동 방정식 환원과 양정부호 시스템#

두 식을 그대로 두면 미지수는 $N_z N_s + N_p$ 개다. 운동량에서 $\mathbf{U}_j^{n+1}$ 를 $\eta^{n+1}$ 의 함수로 풀어 자유표면에 넣는다. $\mathbf{A}$ 가 양정부호면 역행렬도 양정부호다.

대입하면 이산 파동 방정식이 떨어진다.

미지수는 $\eta^{n+1}$ 하나. 크기는 셀 수 $N_p$. 행렬은 대칭이고, 대각이 강하고, 양정부호다. PCG가 맞다. 논문은 큰 격자에서도 반복이 수십 번에서 멈춘다고 적는다.

연직 구조가 없는 1D에서는 $\mathbf{A}$ 가 1이라, 코드의 $k = g H \Delta t^2 \theta^2 / \Delta x^2$ 가 이웃 계수다. $\Delta t$ 가 명시 한계의 4배이고 $c^2 = gH$ 이면 $k = 16\theta^2$ 이다. θ=0.6이면 $k=5.76$. 내부 대각 $1+2k$ 는 12.52, 비대각 둘의 합은 11.52. 대각이 이긴 폭은 언제나 1이고, 그 1은 저장항 $P_i \eta$ 가 남긴 것이다. 스텝을 키우면 $k$ 가 대각과 비대각에 같이 붙어 잉여는 늘지 않는다. 앞 슬라이더의 CFL=6, θ=0.6이면 $k=12.96$, 대각 26.92, 비대각 합 25.92. 폭은 여전히 1이다. 행렬은 더 음해로 가도 더 안전해지지 않는다. 죽는 것은 진폭이다.

$\eta$ 가 풀리면 $\mathbf{U}_j^{n+1}$ 는 변마다 작은 삼중대각으로 떨어진다. 연직 속도는 비압축 식으로 거슬러 올라간다.

Python으로 1D 반음해 풀이#

같은 분리다. $N=200$, $L=1$, $H=1$, $g=9.81$, θ=0.6.

import numpy as np
from scipy.linalg import solve_banded
 
class CasulliShallowWater1D:
    """1D 선형 천수 — 반음해 θ 방법."""
 
    def __init__(self, N=200, L=1.0, H=1.0, g=9.81, theta=0.6):
        self.N, self.L, self.H, self.g, self.theta = N, L, H, g, theta
        self.dx = L / N
        self.eta = np.zeros(N)            # 셀 중심 자유표면
        self.u = np.zeros(N + 1)          # 변(face) 수직속도, 양 끝 = 0
 
    def initial_bump(self, x0=0.3, sigma=0.05, amp=0.05):
        x = (np.arange(self.N) + 0.5) * self.dx
        self.eta = amp * np.exp(-((x - x0) ** 2) / (2 * sigma ** 2))
        self.u[:] = 0.0
 
    def step(self, dt):
        N, dx, g, H, th = self.N, self.dx, self.g, self.H, self.theta
        k = g * H * dt * dt * th * th / (dx * dx)   # 파동 항 계수
        alpha = H * dt / dx
 
        # G_j: 명시적 부분 모은 변값(양 끝 = 0)
        G = np.zeros(N + 1)
        G[1:N] = self.u[1:N] - g * (dt / dx) * (1 - th) * (self.eta[1:] - self.eta[:-1])
 
        # 삼중대각: -k η_{i-1} + (1+2k) η_i - k η_{i+1} = R_i
        ab = np.zeros((3, N))
        ab[0, 1:] = -k                         # super-diagonal
        ab[2, :-1] = -k                        # sub-diagonal
        ab[1, :] = 1.0 + 2 * k
        ab[1, 0] = 1.0 + k                     # 좌벽
        ab[1, -1] = 1.0 + k                    # 우벽
 
        rhs = self.eta.copy()
        rhs -= alpha * th * (G[1:] - G[:-1])
        rhs -= alpha * (1 - th) * (self.u[1:] - self.u[:-1])
 
        new_eta = solve_banded((1, 1), ab, rhs)
 
        new_u = self.u.copy()
        new_u[1:N] = self.u[1:N] - g * (dt / dx) * (
            th * (new_eta[1:] - new_eta[:-1]) + (1 - th) * (self.eta[1:] - self.eta[:-1])
        )
        new_u[0] = new_u[-1] = 0.0
        self.eta, self.u = new_eta, new_u
        return self.eta
 
if __name__ == "__main__":
    sim = CasulliShallowWater1D(N=200, theta=0.6)
    sim.initial_bump()
    c = np.sqrt(sim.g * sim.H)
    dt_explicit_limit = sim.dx / c
    dt = 4.0 * dt_explicit_limit                # 명시 한계의 4배
    for _ in range(int(0.6 / dt)):
        sim.step(dt)
    print(f"max|η| = {np.max(np.abs(sim.eta)):.4f}  (안정 유지)")

명시는 $\Delta t = 4 \Delta t_\text{ex}$ 에서 곧 발산한다. Casulli는 같은 간격으로 구간 0.6을 끝낸다. 위상은 늦고, 진폭은 θ가 1에 가까울수록 더 죽는다. θ=0.5는 가중치가 반반이라 2차다. 0.6을 넘기면 1차다. 코드 기본값 0.6은 그 1차 쪽에 발을 걸친 자리다.

시간 스텝을 키워보자#

수로 길이 1, 셀 200이면 $\Delta x = 0.005$. $c=\sqrt{9.81}\approx 3.13$ m/s라 명시 한계는 약 0.0016초, 네 배는 약 0.0064초다. 0.6을 이 간격으로 나눈 뒤 소수점을 버리면 93스텝이다. 출력은 max|η| 가 유한하다는 한 줄뿐이라, 초기 융기 0.05가 얼마나 줄었는지는 그 줄이 말하지 않는다.

그 0.05는 수심 $H=1$ 의 5%다. 선형화 $\partial_t \eta + H\partial_x u = 0$ 은 $\eta$ 가 $H$ 보다 충분히 작아야 닫힌다. 그래서 셀러리티는 $\sqrt{gH}$ 한 값이고, 삼중대각의 $k$ 도 셀마다 같다. 융기를 수심만큼 키우면 $c$ 가 자리에 따라 달라지고, 이 일정한 $k$ 로는 파동 방정식이 닫히지 않는다.

θ=0.5에 같은 4배 스텝을 넣으면 $k=16\times 0.25=4$ 다. 대각 우위는 역시 1이다. 안정 경계라는 말은 행렬이 아니라 시간 가중 $(θ-1/2)$ 에서 온다. 그 값이 0이면 수치 점성이 꺼지고, 격자 모드가 감쇠 없이 왕복한다.

비판적 고찰#

만능은 아니다.

시간 차수부터 값이 다르다. θ=0.5는 2차지만 안정 경계라 중력파 파면이 쉽게 휜다. 0.6을 넘으면 1차다. 조위 진폭의 숫자가 산물이면 스텝을 줄이거나, 더 높은 차수의 IMEX로 옮긴다.

이송은 명시 항에 남는다. 단순 업윈드는 수치 점성이 세다. Lagrange 궤적을 거꾸로 밟아 보간하는 EL을 붙이는 이유가 그것이다. 격자 한 칸을 넘는 $\Delta t$ 에서도 안정하지만, 보간 차수와 출발점 오차에 결과가 흔들린다.

정수압 $\partial p/\partial z = -\rho g$ 도 전제다. 쇄파나 낙차처럼 연직 가속도가 큰 곳은 비정수압 보정이 필요하다.

둔각 하나면 5분 스텝을 거둔다#

5분과 θ=0.6의 거래는 스텝 수를 1500분의 1로 줄이고, 파고의 1차 감쇠를 받아들이는 일이다. 검조 기록의 진폭이 산물이면 θ를 1 쪽으로 밀지 않는다. 대각이 이긴 폭은 이미 1이라, 더 음해로 가도 행렬이 더 안전해지지 않는다. 죽는 것은 파다.

이 차분은 중심을 잇는 선이 변에 수직일 때만 $\partial \eta / \partial n$ 이다. 둔각 하나가 외심을 셀 밖으로 밀어내면 등식은 끝난다. 면에 수직이 아닌 성분이 남는데, 그 성분을 확산 플럭스와 같은 꼴로 되돌릴 수 있는지는 비직교 면의 확산 보정이 묻는 질문이다.