N체 중력 샌드박스베타

N체 중력 샌드박스 소개

이 샌드박스는 화면에 있는 모든 천체 쌍 사이의 뉴턴 중력을 하나도 빼지 않고 더한다(직접합, 계산량은 천체 수의 제곱에 비례). 그렇게 구한 가속도로 위치와 속도를 조금씩 전진시키는데, 그 방법은 「속도 베를레(kick-drift-kick)」다. 속도를 반 스텝 밀고(kick), 위치를 한 스텝 옮기고(drift), 새 위치에서 힘을 다시 구해 속도를 나머지 반 스텝 민다.

화면 오른쪽의 「에너지 오차」는 시작할 때의 총에너지(운동 + 위치 에너지)와 지금 값의 차이를 백분율로 보여 준다. 중력만 작용하는 계에서 총에너지는 보존되어야 하므로, 이 숫자는 곧 수치 적분의 오차다. 병합이나 천체 추가가 일어나면 기준을 다시 잡는다.

시간 단위와 질량 단위는 게임용으로 정한 값이다(태양급 항성 질량 4000, 중력 상수 기본 100). 궤도의 모양과 보존 법칙은 실제와 같지만, 「초」나 「질량」을 실제 천문 단위로 바꿔 읽을 수는 없다.

왜 오일러 방법으로는 행성이 달아나는가

가장 단순한 적분법인 오일러 방법은 지금의 속도로 위치를, 지금의 가속도로 속도를 한꺼번에 갱신한다. 타원 궤도에서 이 방법은 매 바퀴 에너지를 조금씩 더한다. 오차가 서로 상쇄되지 않고 한 방향으로 쌓이기 때문이다.

직접 계산해 보면 차이가 선명하다. 중심 질량 GM = 1, 거리 1에서 접선 속도 0.8로 출발한 궤도(이심률 0.36, 주기 약 3.96)를 시간 간격 0.01로 적분했다. 오일러 방법은 한 바퀴 뒤 에너지가 22.9% 늘었고, 10바퀴 뒤 57%, 100바퀴 뒤 89%였다. 궤도가 바퀴마다 부풀어 결국 풀려나는 쪽으로 간다. 같은 조건에서 kick-drift-kick은 어느 순간에도 오차가 0.021%를 넘지 않았고, 100바퀴가 지나도 커지지 않았다.

이 차이는 스텝을 줄여서 메울 수 있는 종류가 아니다. 베를레 계열 방법은 「심플렉틱」 적분법이어서, 정확한 에너지 대신 그와 아주 가까운 수정된 에너지를 보존한다. 그래서 오차가 궤도 한 바퀴 안에서 오르내릴 뿐 누적되지 않는다(Hairer·Lubich·Wanner 2003). 분자동역학에서 Verlet(1967)이 이 방법을 쓴 이유도 수십만 스텝 동안 에너지가 새지 않아야 했기 때문이다. 이 도구의 「적분기」를 「오일러 (비교용)」로 바꾸면 같은 궤도에서 「에너지 오차」가 계속 커지고 궤도가 바깥으로 풀려 나가는 것을 직접 볼 수 있다.

에너지 오차 숫자를 키우는 것들

베를레 방법도 스텝이 거칠면 오차가 커진다. 이 도구는 한 프레임(1/60초)을 여러 부분 스텝으로 나누는데, 천체가 적을수록 많이(최대 16) 나눈다. 「시간 배속」을 올리면 한 부분 스텝이 담당하는 시간이 길어져 오차가 커진다. 태양계 프리셋에서 배속을 1과 4로 바꿔 가며 오차 표시를 비교해 보면 된다.

가장 큰 오차는 두 천체가 아주 가까이 스칠 때 생긴다. 힘이 거리의 제곱에 반비례하므로 근접 통과 순간 가속도가 한 스텝 사이에 급변한다. 「소프트닝」은 거리 계산에 작은 길이 ε를 더해(√(r² + ε²)) 이 발산을 누그러뜨린다. 대신 ε보다 가까운 거리에서는 더 이상 뉴턴 중력이 아니다. 소프트닝을 0으로 두고 행성을 항성 바로 옆으로 던지면 비정상적으로 빠른 튕겨 나감이 나오기도 하는데, 이것은 물리가 아니라 수치 오차다.

천체가 화면 밖 아주 먼 곳(좌표 5만)으로 달아나면 제거된다. 병합·파괴·추가·삭제 때와 마찬가지로 천체 수가 바뀌는 순간 기준 에너지를 다시 잡으므로, 오차 표시는 남은 천체들의 에너지를 기준으로 이어진다.

세 천체 문제와 8자 궤도

두 천체의 운동은 케플러 궤도로 정확히 풀리지만, 세 번째가 들어오면 일반적인 닫힌 해가 없다. 대부분의 초기 조건에서 운동은 혼돈적이어서, 처음 위치가 아주 조금만 달라도 얼마 뒤 궤적이 완전히 달라진다. 같은 등급의 항성 세 개를 비슷한 속도로 던져 보면 한동안 얽혀 돌다가 하나가 튕겨 나가고 둘이 쌍성으로 남는 경우가 흔하다.

그런데도 아름다운 주기해가 있다. 질량이 같은 세 천체가 하나의 8자 곡선 위를 서로 쫓아 도는 해다. Moore(1993)가 수치 계산으로 먼저 찾았고, Chenciner와 Montgomery(2000)가 변분법으로 존재를 증명했다. G = m = 1 단위에서 널리 쓰이는 초기 조건(Simó 계산값)은 두 천체가 (±0.97000436, ∓0.24308753), 가운데 천체의 속도가 (−0.93240737, −0.86473146), 바깥 두 천체의 속도는 각각 그 절반의 반대 방향이며 주기는 약 6.3259다.

이 초기 조건을 시간 간격 0.001로 적분해 보면 차이가 다시 드러난다. kick-drift-kick은 10주기 뒤에도 가운데 천체가 원점으로 돌아오고 에너지 오차는 10⁻¹¹% 수준이었다. 오일러 방법은 한 주기 만에 에너지가 3% 어긋났고 10주기 뒤에는 20% 어긋나 8자 모양이 무너졌다. 「8자 궤도 (삼체)」 프리셋은 이 초기 조건을 이 도구의 단위에 맞게 늘려 놓는다. 세 별의 질량 M = 4000, 길이 단위 L은 보이는 영역 짧은 변 절반의 0.42배(화면이 작으면 가장 가까이 올 때도 두 별이 닿지 않을 만큼 더 크게)이고, 속도는 √(G·M/L)배, 주기는 6.3259 × L ÷ √(G·M/L)가 된다. 적분기를 오일러로 바꾸고 같은 프리셋을 다시 불러오면 8자가 무너지는 과정도 비교할 수 있다.

충돌 모델은 근사다

「현실적」 충돌은 충격 에너지와 결합 에너지의 비로 병합, 스치고 지나감, 부분 파괴, 완전 파괴를 나눈다. 질량·운동량·질량 중심은 보존하지만, 파편 수와 크기, 초신성·킬로노바가 일어나는 질량 문턱은 게임용으로 정한 값이다. 행성 형성이나 거대 충돌 연구에서 쓰는 SPH 같은 유체 계산과는 다르다.

충돌 판정은 부분 스텝 동안 두 천체가 곧게 움직였다고 보고, 그 구간에서 가장 가까웠던 거리가 두 반지름의 합보다 작았는지를 본다(선분 최근접 검사). 그래서 아주 빠른 작은 천체도 한 스텝에 상대를 건너뛰어 지나가지 않는다. 서로 지나친 뒤에 잡힌 충돌은 질량 중심은 그대로 두고 가장 가까웠던 순간의 배치로 되돌려 충돌 방향을 정한다.

예시

사용 방법

1. 프리셋(태양계, 쌍성계, 원시 원반 등)을 고르거나 빈 화면에서 시작한다.
2. 천체 종류를 고른 뒤 끌었다 놓아 던진다. 톡 누르면 가장 무거운 천체를 도는 원궤도 속도로 놓인다.
3. 중력 상수·시간 배속·소프트닝·적분기(속도 베를레/오일러)를 바꾸면서 궤도와 「에너지 오차」가 어떻게 변하는지 본다.
4. 필요하면 3D로 바꿔 궤도면이 기울어진 계를 돌려 본다.

주요 용어

심플렉틱 적분법

위상 공간의 부피를 보존하는 적분법. 에너지 오차가 쌓이지 않고 오르내린다. 속도 베를레, 립프로그가 대표적이다.

소프트닝(ε)

중력 계산에서 거리 r 대신 √(r² + ε²)을 쓰는 것. 근접 통과 때 힘이 무한대로 커지는 것을 막는다.

안무 해(choreography)

여러 천체가 같은 닫힌 곡선을 시간 차를 두고 따라가는 주기해. 8자 궤도가 첫 사례로 유명하다.

자주 묻는 질문

참고

2D 모드에서는 모든 천체가 한 평면(z = 0)에 있고 z 방향 힘이 0이 된다. 3D에서 2D로 바꾸면 천체가 평면으로 눌리므로 궤도가 달라진다.