SPH

1. 개요

SPH Smoothed Particle Hydrodynamics
정식 명칭	평활 입자 유체동역학
분야	전산유체역학 × 컴퓨터 그래픽스 × 천체물리학
기법 분류	무격자(mesh-free) · 라그랑주(Lagrangian)
기원	1977년, 천체물리 시뮬레이션
핵심 개념	커널 함수, 평활 길이
대표 커널	3차 스플라인(지지 반경 2h) · 그래픽스계 poly6/spiky(지지 반경 h)
이웃 수	3차원에서 h = 1.2Δx 이면 56개(연속체 추정 57.9)

격자를 짜기 싫으면, 격자를 안 쓰면 된다. 대신 유체를 수십만 개의 점으로 부순다.

SPH(Smoothed Particle Hydrodynamics, 평활 입자 유체동역학)는 유체를 고정된 격자 대신 움직이는 입자(particle)의 집합으로 표현하는 무격자(mesh-free) 라그랑주(Lagrangian) 수치 기법이다. 각 입자는 질량·속도·밀도 같은 물성을 지닌 채 유체를 따라 흘러가며, 임의 위치의 물리량은 주변 입자들의 값을 커널 함수로 가중 평균해 추정한다. 격자 생성이라는 전산유체역학 최대의 고행에서 해방된다는 것이 이 기법의 가장 큰 매력이다.

물이 흩뿌려지고, 파도가 부서지고, 유체가 갈라졌다 다시 합쳐지는 것처럼 형상이 극적으로 변하고 자유표면이 요동치는 문제에서 SPH는 빛을 발한다. 격자 기반 기법은 이런 문제에서 격자를 매 순간 다시 짜거나 계면을 추적하느라 고생하지만, 입자는 그냥 유체를 따라 흘러가면 그만이기 때문. 영화 속 쏟아지는 물, 게임 속 출렁이는 용암은 상당수가 이 계열의 입자법으로 태어난다.¹

유체역학의 지배식을 입자로 옮긴다는 점에서 SPH는 물질점법·FLIP-PIC 유체 같은 입자 기반 기법들과 한 가족이고, 컴퓨터 그래픽스 쪽에서는 사실상 이 계열의 입문 관문 역할을 한다.

2. 역사: 별에서 시작한 물

의외로 SPH의 고향은 유체공학이 아니라 천체물리학이다. 1977년, 별과 은하처럼 텅 빈 우주 공간을 떠도는 물질 덩어리를 시뮬레이션하려는 문제에서 태어났다. 이런 문제에는 미리 격자를 깔 “공간”이라는 것 자체가 무의미했고, 그래서 물질 자체를 입자로 보는 발상이 자연스러웠다.

• 1977년 — 진 몬해건(Joe Monaghan)과 밥 깅골드(Bob Gingold)가, 그리고 별개로 리언 루시(Leon Lucy)가 거의 동시에 SPH를 제안한다.² 각각 은하의 진동(비구형 별의 안정성)과 별의 분열 가설 검증을 다루고 있었다.
• 1980~90년대 — 초신성 폭발, 은하 충돌, 행성 형성 같은 천체물리 시뮬레이션의 주력 도구로 자리 잡는다.
• 1992년 — 몬해건의 Annual Review 리뷰가 흩어져 있던 공식들을 한 권으로 정리한다. 커널·인공점성·시간 간격의 “국룰”이 사실상 이 논문에서 굳었다.
• 1994년 — 몬해건이 자유표면 유동에 SPH를 적용하며 지상의 유체공학으로 무대를 넓힌다. 오늘날 약압축성 SPH(WCSPH)라 불리는 노선의 출발점이다.
• 2003년 — 뮐러(Müller) 등이 실시간 인터랙티브 유체를 위한 SPH를 컴퓨터 그래픽스 학계에 소개하며, 게임·영화 VFX로 급속히 퍼진다.
• 2009~2014년 — PCISPH·IISPH처럼 비압축성을 압력 반복으로 강제하는 개량판이 그래픽스 쪽에서 쏟아진다.

별을 흩뿌리던 방법이 물을 흩뿌리는 방법이 된 셈이다.

3. 커널과 평활 길이

SPH의 수학적 심장은 커널 함수(smoothing kernel) $W$다. 어떤 물리량 $A$의 위치 $\mathbf{r}$에서의 값은, 주변 입자 $j$들의 값을 커널로 가중해 다음처럼 근사한다.

$$A(\mathbf{r}) \approx \sum_{j} m_j \frac{A_j}{\rho_j} \, W(\mathbf{r} - \mathbf{r}_j, h)$$

여기서 $m_j$는 입자 $j$의 질량, $\rho_j$는 밀도, $h$는 평활 길이(smoothing length)다. 평활 길이가 “이웃”으로 칠 범위를 정한다 — $h$가 작으면 국소적이지만 이웃이 적어 잡음이 커지고, $h$가 크면 매끄럽지만 계산이 무거워지고 세부가 뭉개진다. 이 절충이 SPH 정확도의 급소다.

커널이 만족해야 할 조건은 네 가지다. 규격화($\int W \, dV = 1$), 대칭성(짝함수 — 이 덕에 근사가 1차가 아니라 2차가 된다), 콤팩트 지지(유한 반경 밖에서 0 — 이게 없으면 이웃 탐색이 $O(N^2)$로 폭발한다), 그리고 $h \to 0$에서 디랙 델타로 수렴.

여기서 입문자를 가장 자주 넘어뜨리는 지점이 지지 반경과 $h$의 관계다. 두 관례가 공존하기 때문이다.

• 물리·천체물리 관례(몬해건 계열) — 3차 스플라인 커널의 지지 반경은 $2h$다. 즉 $h$는 “잘리는 거리”가 아니라 표준편차에 가까운 폭 척도이고, 실제로 이웃을 끊는 거리는 그 두 배다.
• 그래픽스 관례(뮐러 2003 계열) — poly6·spiky 커널은 $r > h$에서 정확히 0이다. 여기서는 $h$가 곧 지지 반경이다.

같은 기호 $h$가 경우에 따라 두 배 차이 나는 것을 모르고 코드를 섞으면, 이웃 수가 8배(3차원) 틀어지면서 밀도가 조용히 어긋난다. 논문을 읽을 때 "$\kappa h$"라는 표기가 보이면 $\kappa$가 바로 이 비례상수다.

4. 커널 규격화 — 직접 검산

가장 널리 쓰이는 3차 스플라인(cubic spline, $M_4$) 커널은 $q = r/h$로 쓰면 이렇다.

$$W(r,h) = \sigma_d \begin{cases} 1 - \tfrac{3}{2}q^2 + \tfrac{3}{4}q^3, & 0 \le q < 1 \\ \tfrac{1}{4}(2-q)^3, & 1 \le q < 2 \\ 0, & q \ge 2 \end{cases}$$

차원별 규격화 상수 $\sigma_d$는 암기 사항처럼 인용되지만, 실은 $\int W \, dV = 1$에서 바로 떨어진다. 직접 적분해 보면 이렇게 맞는다.

차원	규격화 상수	필요한 무차원 적분	적분값(수치)
1D	2/(3h)	∫ f dq (0→2) = 3/4	0.750000
2D	10/(7πh²)	∫ q f dq (0→2) = 7/20	0.350000
3D	1/(πh³)	∫ q² f dq (0→2) = 1/4	0.250000

가령 3차원에서는 $\int W \, dV = \sigma_3 \cdot 4\pi h^3 \int_0^2 q^2 f(q)\,dq = \frac{1}{\pi h^3}\cdot 4\pi h^3 \cdot \frac14 = 1$로 정확히 닫힌다. 1D·2D도 같은 식으로 딱 맞는다(심프슨 적분 200만 점으로 세 값 모두 소수점 12자리까지 재현). 즉 $10/(7\pi)$ 같은 괴상한 상수는 누가 외우라고 만든 게 아니라 적분이 토해낸 숫자다.

같은 커널의 2차 모멘트도 닫힌 값을 준다. 3차원에서

$$\langle r^2 \rangle = \int r^2 W \, dV = 0.9\,h^2 \quad \Rightarrow \quad r_\text{rms} = 0.9487\,h$$

이 숫자가 왜 중요한가 하면, 유효 해상도는 $h$가 아니라 $r_\text{rms}$ 쪽에 가깝다는 뜻이기 때문이다. $h$만 보고 “격자 간격 $h$짜리 해상도”라고 말하면 과장이 된다.

5. poly6·spiky — 그래픽스 커널의 함정

실시간 그래픽스에서 애용되는 뮐러 2003의 커널 두 개도 같은 방식으로 검산된다.

$$W_\text{poly6} = \frac{315}{64\pi h^9}(h^2-r^2)^3, \qquad W_\text{spiky} = \frac{15}{\pi h^6}(h-r)^3 \qquad (0 \le r \le h)$$

두 상수 모두 수치적분으로 확인된다($1/(4\pi\int_0^1 x^2(1-x^2)^3 dx) = 1.566681 = 315/64\pi$, $1/(4\pi\int_0^1 x^2(1-x)^3 dx) = 4.774648 = 15/\pi$). poly6의 미덕은 $r^2$만 쓰므로 제곱근을 한 번도 안 뽑는다는 것 — 실시간에서는 이게 곧 성능이다.

그런데 poly6에는 유명한 함정이 있다. 기울기가

$$\nabla W_\text{poly6} \propto -r\,(h^2-r^2)^2$$

라서 $r \to 0$에서 0으로 사라진다. 입자가 겹치다시피 몰려도 밀어내는 힘이 안 생긴다는 뜻이다. $h=1$로 두고 실제로 재면 $r = 0.01h$에서 $|\partial W/\partial r|$이 poly6는 0.094, spiky는 14.04다 — 150배 차이. 그래서 뮐러 등은 밀도에는 poly6를, 압력 기울기에는 기울기가 죽지 않는 spiky를 쓰고, 점성에는 라플라시안이 음수로 발산하지 않는 또 다른 커널을 쓴다. “커널 하나로 통일”이 아니라 항마다 다른 커널이라는 그래픽스 SPH의 희한한 관례가 여기서 나왔다.

6. 밀도·압력·힘의 근사

SPH가 아름다운 이유는, 미분 연산을 커널의 미분으로 떠넘길 수 있다는 데 있다. 즉 유체의 기울기나 발산을 구할 때 물리량을 미분하는 대신 커널 함수를 미분하면 된다 — 그리고 커널은 우리가 고른 매끄러운 함수라 미분이 해석적으로 알려져 있다.

이 원리로 나비에-스토크스 방정식의 각 항이 입자합으로 번역된다.

1. 밀도 — 각 입자 주변의 이웃 질량을 커널로 합산해 $\rho_i = \sum_j m_j W_{ij}$로 구한다.
2. 압력 — 밀도로부터 상태방정식으로 계산한다. 밀도가 기준보다 높으면(입자가 몰리면) 압력이 커져 서로를 밀어낸다.
3. 힘 — 압력 기울기 항과 점성 항을 각각 커널의 기울기·라플라시안으로 근사해 입자에 작용하는 가속도를 얻는다.

이렇게 얻은 가속도로 입자를 시간적분하면(대개 베를레 적분이나 심플렉틱 적분기) 다음 스텝의 위치와 속도가 나온다. 격자도, 연립방정식도 없이 오로지 이웃 탐색과 커널 합만으로 유체가 흐른다. 실무 비용의 절반 이상이 그 이웃 탐색에 들어가므로, 균일 격자 해시나 최근접 이웃 탐색 자료구조가 사실상 필수다.

점성 항은 실제 점성 계수를 그대로 쓰기보다 인공점성(artificial viscosity)을 쓰는 쪽이 천체물리에서 표준이다. 충격파를 소산으로 잡아 주는 장치인데, 그 대가로 매끈한 영역에도 소산이 깔린다. 3차원 3차 스플라인에서 인공점성 계수 $\alpha$는 대략

$$\nu_\text{eff} \approx 0.1\,\alpha\,c\,h$$

정도의 실제 운동점성에 해당한다.³ 이 식의 함의가 아프다 — 달성 가능한 레이놀즈수가 해상도에 묶인다. 특성길이 $L$, 속도 $V$, $c = 10V$, $h = L/100$, $\alpha = 0.01$을 넣으면 $\nu_\text{eff} = 10^{-4} VL$이므로 $Re \approx 10^4$가 상한이다. 더 높은 $Re$를 보고 싶으면 입자를 더 잘게 쪼개는 길밖에 없다.

7. 이웃 수와 격자 밀도 검산

실무에서 $h$는 입자 간격 $\Delta x$에 비례해 잡는다($h = \eta \Delta x$, 보통 $\eta = 1.2 \sim 1.5$). 지지 반경이 $2h$이므로 이웃 수는 기하학적으로 정해진다.

$$N_\text{neigh} \approx \frac{4}{3}\pi (2\eta)^3 \quad (\text{3D}), \qquad \pi (2\eta)^2 \quad (\text{2D})$$

$\eta = 1.2$면 3차원 57.9, 2차원 18.1이다. 정육면체 격자에 실제로 입자를 깔고 세어 보면 56개(2차원 20개)로, 연속체 추정과 몇 개 차이로 맞는다.

더 재미있는 것은 같은 배치에서 밀도합이 얼마나 정확한가다. 완벽한 격자에 $m = \rho_0 \Delta x^3$인 입자를 깔고 $\rho_i = \sum_j m_j W_{ij}$를 그대로 돌리면 이렇게 나온다.

h/Δx	이웃 수(3D 격자)	밀도합 ρ/ρ₀	오차
1.0	26	0.999972	−0.003%
1.2	56	1.000810	+0.081%
1.3	80	0.997262	−0.274%
1.5	92	1.001795	+0.180%

정확히 1이 되지 않는다. 그리고 더 중요한 것은 오차가 $\eta$에 따라 단조롭게 줄지 않는다는 점이다 — $\eta=1.3$이 $\eta=1.2$보다 이웃을 24개 더 쓰면서 오차는 3배 크다. 커널 합이 격자점과 어떻게 맞물리는지(격자 공명)에 따라 부호까지 바뀌기 때문이다. 그래서 SPH에서 “이웃을 늘리면 정확해진다”는 직관은 반만 맞고, 코드마다 $\eta$가 1.2, 1.3, 1.5로 갈리는 이유도 여기에 있다. 밀도의 절대값이 중요한 문제라면 초기 배치에서 $\rho/\rho_0$를 직접 재 보고 질량을 보정하는 것이 안전하다.

8. 약압축성(WCSPH)과 시간 간격

물처럼 거의 안 눌리는 유체를 SPH로 다루는 가장 흔한 길은 약압축성(weakly compressible) 접근이다. 테이트(Tait) 상태방정식

$$p = B\left[\left(\frac{\rho}{\rho_0}\right)^{\gamma} - 1\right], \qquad \gamma = 7, \quad B = \frac{\rho_0 c^2}{\gamma}$$

으로 압력을 밀도에 직결시킨다. 여기서 $c$는 실제 음속이 아니라 우리가 고른 수치적 음속이다. 선형화하면 $p \approx c^2 \delta\rho$이고, 동압 $\rho_0 V^2/2$를 떠받치려면

$$\frac{\delta\rho}{\rho_0} \approx \frac{V^2}{2c^2} = \frac{M^2}{2}$$

만큼 밀도가 출렁여야 한다(마하수 $M = V/c$). 몬해건의 국룰인 $c \ge 10 V$가 바로 여기서 나온다 — $M = 0.1$이면 밀도 요동이 0.5%로 억제된다.

왜 실제 음속을 안 쓰느냐는 CFL 조건이 답한다. 음속 기반 시간 간격은 $\Delta t \lesssim 0.25\,h/c$ 꼴이고, 물의 실제 음속 1480 m/s와 $h = 1$ cm를 넣으면

음속 c	B = ρ₀c²/γ	마하수(V=2 m/s)	밀도 요동	Δt 한도
1480 m/s (실제)	3.13×10⁸ Pa	0.00135	0.0001%	1.69×10⁻⁶ s
20 m/s (c=10V)	5.71×10⁴ Pa	0.1	0.5%	1.25×10⁻⁴ s
10 m/s (c=5V)	1.43×10⁴ Pa	0.2	2%	2.50×10⁻⁴ s

74배다. 물을 “조금 눌리게” 타협하는 대가로 스텝 수를 두 자릿수 줄이는 셈이고, 이것이 WCSPH가 살아남은 이유 전부다.⁴ 반대로 비압축성을 진짜로 강제하려면 압력 푸아송 방정식을 풀거나(ISPH), 밀도 오차를 반복으로 눌러야 한다(PCISPH·IISPH) — 스텝당 비용이 오르는 대신 $\Delta t$를 훨씬 키운다.

9. 격자법과의 비교, 그리고 그 대가

SPH의 장점과 단점은 동전의 양면이다.

항목	SPH (입자·라그랑주)	격자법 (FVM/FDM·오일러)
격자 생성	불필요	필수, 종종 최대 난관
자유표면·큰 변형	자연스럽게 처리	계면 추적 필요
질량 보존	입자 수로 정확히 보존	이산화에 의존
비압축성	까다로움	성숙한 기법 다수
경계 조건	골칫거리	잘 정립됨
정확도(동일 비용)	대체로 낮음	대체로 높음

가장 아픈 급소는 두 가지다. 첫째, 비압축성 — 위에서 본 음속·시간 간격 절충이 그것이다. 둘째, 경계 조건 처리. 벽 근처에서는 이웃 입자가 한쪽에만 있어 커널 합이 편향된다. 위 격자 검산에서 $\rho/\rho_0 \approx 1$이 나온 것은 입자가 사방에 있었기 때문이고, 벽에서는 지지구의 절반이 비어 밀도가 곧바로 절반 가까이 내려앉는다(입자 결핍, particle deficiency). 이를 메우려 유령 입자·경계 입자·커널 보정 같은 온갖 땜질이 동원된다.

게임·영화에서는 SPH만 쓰기보다, 입자의 장점과 격자의 장점을 섞은 FLIP-PIC 유체 계열 하이브리드가 사실상 표준으로 자리 잡았다. 입자로 이류(advection)를 처리하고 격자로 압력 투영을 풀어, SPH의 비압축성 약점을 격자로 메우는 전략이다. 상용 툴 대부분이 이 노선을 따른다.

10. 관련 문서

• 나비에-스토크스 방정식 · 유체역학 · 전산유체역학
• 다상유동 · 유한체적법 · 격자
• FLIP-PIC 유체 · 물질점법 · 물리 엔진
• 베를레 적분 · 심플렉틱 적분기 · CFL 조건
• 최근접 이웃 탐색 · 경계 조건 · 마하수 · 점성

11. Footnotes

1. 다만 “물처럼 보이는 것”과 “물리적으로 정확한 물”은 다르다. VFX 아티스트는 대개 정확성보다 화면빨을 택하며, 이 점에서 SPH는 물리 엔진과 같은 실용주의 진영에 속한다. ↩

2. Gingold, R. A. & Monaghan, J. J. (1977), MNRAS 181, 375; Lucy, L. B. (1977), AJ 82, 1013. 같은 해 같은 아이디어가 독립적으로 두 번 발표된, 과학사의 흔한 동시 발견 사례다. 누가 원조냐는 논쟁은 여전히 학회 뒤풀이의 안주다. 참고로 “SPH”라는 약칭을 제목에 박은 쪽은 깅골드-몬해건이고, 루시 논문의 제목에는 그 말이 없다 — 작명이 인용 수를 만든다. ↩

3. 계수 0.1은 3차원 3차 스플라인 + 표준 몬해건 인공점성에 대한 값이고, 커널·차원·인공점성 형태(신호속도형 등)에 따라 대략 1/20~1/10 사이로 갈린다. 그러니 이 식은 “자릿수를 잡는 용도”로만 쓰고, 정량적 점성이 필요하면 물리 점성항을 명시적으로 넣는 쪽이 맞다. ↩

4. 약압축성에서 음속을 낮춰 잡는 것은 “물은 원래 안 눌린다”는 상식을 버리는 일이다. 대신 밀도가 몇 % 출렁이는 것을 받아들인다. 물은 원래 안 눌린다고 배웠지만, 시뮬레이션 앞에서는 다들 조금씩 눌린다. ↩