태양계 컴퓨터시뮬레이션학과 2016년 봄학기 담당교수 : 이형원 E304호, hwlee@inje.ac.kr 운동시뮬레이션 제6주 태양계 컴퓨터시뮬레이션학과 2016년 봄학기 담당교수 : 이형원 E304호, hwlee@inje.ac.kr

다음주 과제 실습해오기 제 5 장 읽어오기

제4장 태양계 케플러(Kepler)의 법칙 역 2제곱 법칙과 행성궤도의 안정성 수성 근일점의 세차 삼체 문제와 목성의 지구에 대한 효과 태양계의 공명 : 커크우드(Kirkwood) 틈과 행성 고리 하이퍼리온(Hyperion)의 카오스 운동

소개 지구상에의 운동은 항상 저항, 공기 저항, 이 존재한다. 행성계는 이상적인 저항 없는 운동이다.

케플러 법칙 𝑟 = 𝑥,𝑦 뉴톤의 중력 법칙 태양의 운동은 무시 운동 방정식 𝐹 𝐺 = 𝐺 𝑀 𝑠 𝑀 𝐸 𝑟 2 𝐹 𝐺 𝐹 𝑥 𝐹 𝑦 𝜃 𝑟 = 𝑥,𝑦 뉴톤의 중력 법칙 𝐹 𝐺 = 𝐺 𝑀 𝑠 𝑀 𝐸 𝑟 2 태양의 운동은 무시 운동 방정식 𝑑 2 𝑥 𝑑 𝑡 2 = 𝐹 𝐺,𝑥 𝑀 𝐸 , 𝑑 2 𝑦 𝑑 𝑡 2 = 𝐹 𝐺,𝑦 𝑀 𝐸 𝐹 𝐺,𝑥 =− 𝐺 𝑀 𝑆 𝑀 𝐸 𝑟 2 cos 𝜃 =− 𝐺 𝑀 𝑆 𝑀 𝐸 𝑥 𝑟 3 𝐹 𝐺,𝑦 =− 𝐺 𝑀 𝑆 𝑀 𝐸 𝑟 2 sin 𝜃 =− 𝐺 𝑀 𝑆 𝑀 𝐸 𝑦 𝑟 3

태양계의 특징 행성 반경(AU) 질량(Kg) 이심률 태양(Sun) … 2.0× 10 30 수성(Mercury) 0.39 2.4× 10 23 0.206 금성(Venus) 0.72 4.9× 10 24 0.007 지구(Earth) 1.00 6.0× 10 24 0.017 화성(Mars) 1.52 6.6× 10 23 0.093 목성(Jupiter) 5.20 1.9× 10 27 0.048 토성(Saturn) 9.54 5.7× 10 26 0.056 천왕성(Uranus) 19.19 8.8× 10 25 0.046 해왕성(Neptune) 30.06 1.03× 10 26 0.010 명왕성(Pluto) 39.53 ~6.0×4 0.248

일차 미분 방정식 𝑑 𝑣 𝑥 𝑑𝑡 =− 𝐺 𝑀 𝑠 𝑥 𝑟 3 𝑑𝑥 𝑑𝑡 = 𝑣 𝑥 𝑑 𝑣 𝑦 𝑑𝑡 =− 𝐺 𝑀 𝑠 𝑦 𝑟 3 𝑑 𝑣 𝑥 𝑑𝑡 =− 𝐺 𝑀 𝑠 𝑥 𝑟 3 𝑑𝑥 𝑑𝑡 = 𝑣 𝑥 𝑑 𝑣 𝑦 𝑑𝑡 =− 𝐺 𝑀 𝑠 𝑦 𝑟 3 𝑑𝑦 𝑑𝑡 = 𝑣 𝑦 𝑟~1.5× 10 11 m, 𝑀 𝑠 ≅1.99× 10 30 kg

천문학적 단위 AU : 태양과 지구사이의 평균 거리 (≈1.5× 10 11 m) 1 year (≈3.2× 10 7 s) 지구는 거의 원운동을 한다. 𝑀 𝐸 𝑣 2 𝑟 = 𝐹 𝐺 = 𝐺 𝑀 𝑠 𝑀 𝐸 𝑟 2 (구심력=만유인력) 𝐺 𝑀 𝑠 = 𝑣 2 𝑟= 2𝜋𝑟 1𝑦𝑒𝑎𝑟 2 𝑟=4 𝜋 2 AU 3 / yr 2

천문학적 단위의 미분 방정식 𝑑 𝑣 𝑥 𝑑𝑡 =− 4 𝜋 2 𝑥 𝑟 3 𝑑𝑥 𝑑𝑡 = 𝑣 𝑥 𝑑 𝑣 𝑥 𝑑𝑡 =− 4 𝜋 2 𝑥 𝑟 3 𝑑𝑥 𝑑𝑡 = 𝑣 𝑥 𝑑 𝑣 𝑦 𝑑𝑡 =− 4 𝜋 2 𝑦 𝑟 3 𝑑𝑦 𝑑𝑡 = 𝑣 𝑦

Modelica 방정식 𝑑𝑒𝑟 𝑣𝑥 =−4∗𝑃 𝐼 2 ∗𝑥/𝑟^3; 𝑑𝑒𝑟 𝑥 =𝑣𝑥; 𝑑𝑒𝑟 𝑣𝑦 =−4∗𝑃 𝐼 2 ∗𝑦/𝑟^3; 𝑑𝑒𝑟 𝑦 =𝑣𝑦;

클래스 작성 Motion.y2015.Week06.Earth

초기조건 설정 원 운동을 하는 조건은 구심력과 만유인력이 같아지는 조건이다. 𝑣 2 𝑟 = 4 𝜋 2 𝑟 2 (천문학적 단위 사용) 𝑟=(1,0), 𝑣=(0,2𝜋) 𝑟 𝑣

결과 보이기 보고자 하는 변수 선택

결과 보이기

새 시뮬레이션 결과 타이틀은 Setup으로 수정

결과 해석 초기속도를 2𝜋로 한 경우에는 𝑥,𝑦의 범위가 동일한 원운동을 한다. 초기속도를 4.5로 한 경우에는 𝑥,𝑦의 범위가 동일하지 않은 타원운동을 한다.

케플러 법칙 모든 행성은 태양을 초점으로 하는 타원궤도 운동을 한다. 행성과 태양을 연결한 직선은 같은 시간에 같은 면적을 지난다. 𝑇가 주기이고, 𝑎가 장단축의 길이라고 하면, 𝑇 2 / 𝑎 3 는 상수 이다. Venus: 0.997, Earth: 0.998, Mars: 1.005, Jupiter: 1.010, Saturn: 0.988

Eccentricity

Circle

Ellipse 𝑥= cos 𝑡 , 𝑦= 1− 𝑒 2 sin 𝑡 𝑥 2 + 𝑦 2 1−𝑒 2 2 =1

Parabola 𝑥= 𝑡 2 4𝑎 , 𝑦=𝑡 𝑦 2 =4𝑎𝑥

Hyperbola 𝑥= cosh 𝑡 , 𝑦= 𝑒 2 −1 sinh 𝑡 𝑥 2 − 𝑦 2 𝑒 2 −1 2 =1

역 2제곱 법칙과 행성궤도의 안정성 해석적인 해 𝑟 = 𝑟 2 − 𝑟 1 , 𝜇= 𝑚 1 𝑚 2 𝑚 1 + 𝑚 2 𝑟 = 𝑟 2 − 𝑟 1 , 𝜇= 𝑚 1 𝑚 2 𝑚 1 + 𝑚 2 𝑑 2 𝑑 𝜃 2 1 𝑟 + 1 𝑟 =− 𝜇 𝑟 2 𝐿 2 𝐹 𝑟 𝐹 𝑟 =− 𝐺 𝑀 𝑠 𝑀 𝐸 𝑟 2 𝑟= 𝐿 2 𝜇𝐺 𝑀 𝑠 𝑀 𝐸 1 1−𝑒 cos 𝜃

초기조건 근일점(Perihelion): 𝑟 𝑚𝑖𝑛 =𝑎 1−𝑒 원일점(Aphelion): 𝑟 𝑚𝑎𝑥 =𝑎 1+𝑒 𝑎= 𝐿 2 𝜇𝐺 𝑀 𝑠 𝑀 𝐸 (1− 𝑒 2 ) , 𝐿= 𝜇𝐺 𝑀 𝑠 𝑀 𝐸 𝑎(1− 𝑒 2 ) 𝐿=𝜇 𝑟 𝑚𝑖𝑛 𝑣 𝑚𝑎𝑥 =𝜇 𝑟 𝑚𝑎𝑥 𝑣 𝑚𝑖𝑛 𝑣 𝑚𝑎𝑥 = 𝐺 𝑀 𝑆 1+𝑒 𝑎 1−𝑒 1+ 𝑀 𝐸 𝑀 𝑆 𝑣 𝑚𝑖𝑛 = 𝐺 𝑀 𝑆 1−𝑒 𝑎 1+𝑒 1+ 𝑀 𝐸 𝑀 𝑆

역 제곱 법칙 역 제곱 법칙은 정확한 것인가? 2가 아니고 다른 값이면 어떻게 되는가? 힘선을 고려하면 2가 정확하다.

2가 아닌 경우의 궤도 미분 방정식 일차 미분 방정식 𝐹 𝐺 = 𝐺 𝑀 𝑠 𝑀 𝐸 𝑟 𝛽 𝐹 𝐺 = 𝐺 𝑀 𝑠 𝑀 𝐸 𝑟 𝛽 일차 미분 방정식 𝑑 𝑣 𝑥 𝑑𝑡 =− 4 𝜋 2 𝑥 𝑟 𝛽+1 𝑑𝑥 𝑑𝑡 = 𝑣 𝑥 𝑑 𝑣 𝑦 𝑑𝑡 =− 4 𝜋 2 𝑦 𝑟 𝛽+1 𝑑𝑦 𝑑𝑡 = 𝑣 𝑦

Modelica 코드

PlanetBeta

PlanetBeta

PlanetBeta

PlanetBeta

결과 해석 이심율이 크고 𝛽가 2가 아닌 경우에는 궤도가 안정적이지 않다.

수성 근일점의 세차

관측과의 차이 세차는 다른 행성의 영향으로 발생한다. 566 arcsec/100yr, 1 arcsec = 1/3600 degree 섭동계산: 523 arcsec/100yr 일반상대론: 43 arcsec/100yr 𝐹 𝐺 ≈ 𝐺 𝑀 𝑠 𝑀 𝐸 𝑟 2 1+ 𝛼 𝑟 2 𝛼≈1.1× 10 −8 AU 2 , 𝑎=0.39 AU 𝑟 1 =𝑎 1+𝑒 =0.47 AU, 𝑒=0.206

미분 방정식 𝑑 𝑣 𝑥 𝑑𝑡 =− 4 𝜋 2 𝑥 𝑟 3 1+ 𝛼 𝑟 2 𝑑𝑥 𝑑𝑡 = 𝑣 𝑥 𝑑 𝑣 𝑥 𝑑𝑡 =− 4 𝜋 2 𝑥 𝑟 3 1+ 𝛼 𝑟 2 𝑑𝑥 𝑑𝑡 = 𝑣 𝑥 𝑑 𝑣 𝑦 𝑑𝑡 =− 4 𝜋 2 𝑦 𝑟 3 1+ 𝛼 𝑟 2 𝑑𝑦 𝑑𝑡 = 𝑣 𝑦

수성의 초기 조건 𝑟 1 =𝑎 1+𝑒 =0.47, 𝑣 1 =2𝜋 1−𝑒 𝑎(1+𝑒) 𝑒=0.206, 𝑎=0.39 𝑟 1 =𝑎 1+𝑒 =0.47, 𝑣 1 =2𝜋 1−𝑒 𝑎(1+𝑒) 𝑒=0.206, 𝑎=0.39 Motion.y2015.Week06.PlanetGR

Modelica 코드

결과 해석 𝛼의 값이 적을 수록 수성의 근일점이 천천히 회전하는 것을 알 수 있다.

장단축의 위치 찾기 장단축의 끝을 지날 때 거리가 줄어 든다 주어진 𝛼에 대하여 장단축이 움직이는 각속도 측정 𝑑𝑟 𝑑𝑡 =0, 𝑑 2 𝑟 𝑑 𝑡 2 <0 을 만족한다. 주어진 𝛼에 대하여 장단축이 움직이는 각속도 측정 최소제곱법으로 기울기 계산

𝛼에 따른 세차 각속도 𝛼에 따른 세차 각속도 측정 외삽을 통하여 𝛼=1.1× 10 −8 일 때의 세차 각속도 추정

𝑑𝑟 𝑑𝑡 , 𝑑 2 𝑟 𝑑 𝑡 2 의 계산 𝑟= 𝑥 2 + 𝑦 2 , 𝑣= 𝑣 𝑥 2 + 𝑣 𝑦 2 𝑑𝑟 𝑑𝑡 , 𝑑 2 𝑟 𝑑 𝑡 2 의 계산 𝑟= 𝑥 2 + 𝑦 2 , 𝑣= 𝑣 𝑥 2 + 𝑣 𝑦 2 𝑑𝑟 𝑑𝑡 = 𝑥 𝑣 𝑥 +𝑦 𝑣 𝑦 𝑟 𝑑 2 𝑟 𝑑 𝑡 2 = 𝑣 2 𝑟 − 4𝜋 𝑟 1+ 𝛼 𝑟 2 − 𝑥 𝑣 𝑥 +𝑦 𝑣 𝑦 2 𝑟 3

삼체 문제와 목성의 지구에 대한 효과 𝐹 𝐸𝐽 𝐹 𝐸𝑆 𝐹 𝐸𝐽 지구 목성 𝐹 𝐸𝑆 태양 𝐹 𝐸𝐽 =− 𝐺 𝑀 𝐸 𝑀 𝐽 𝑟 𝐸𝐽 3 𝑟 𝐸𝐽 , 𝑟 𝐸𝐽 = 𝑟 𝐸 − 𝑟 𝐽 𝐹 𝐸 = 𝐹 𝐸𝐽 + 𝐹 𝐸𝑆 , 𝐹 𝐽 = 𝐹 𝐽𝐸 + 𝐹 𝐽𝑆

운동 방정식 𝑀 𝐸 𝑑 2 𝑟 𝐸 𝑑 𝑡 2 = 𝐹 𝐸 , 𝑀 𝐽 𝑑 2 𝑟 𝐽 𝑑 𝑡 2 = 𝐹 𝐽 𝑀 𝐸 𝑑 2 𝑟 𝐸 𝑑 𝑡 2 = 𝐹 𝐸 , 𝑀 𝐽 𝑑 2 𝑟 𝐽 𝑑 𝑡 2 = 𝐹 𝐽 𝐺 𝑀 𝐽 =𝐺 𝑀 𝑆 𝑀 𝐽 𝑀 𝑆 =4𝜋 𝑀 𝐽 𝑀 𝑆 (천문학적 단위) 𝑑 𝑟 𝐸 𝑑𝑡 = 𝑣 𝐸 , 𝑑 𝑣 𝐸 𝑑𝑡 = 𝐹 𝐸 𝑀 𝐸 𝑑 𝑟 𝐽 𝑑𝑡 = 𝑣 𝐽 , 𝑑 𝑣 𝐽 𝑑𝑡 = 𝐹 𝐽 𝑀 𝐽

성분별 운동 방정식 𝑑 𝑥 𝐸 𝑑𝑡 = 𝑣 𝐸𝑥 , 𝑑 𝑣 𝐸𝑥 𝑑𝑡 =− 𝐺 𝑀 𝑆 𝑥 𝐸 𝑟 𝐸 3 − 𝐺 𝑀 𝐽 𝑥 𝐸 − 𝑥 𝐽 𝑟 𝐸𝐽 3 𝑑 𝑦 𝐸 𝑑𝑡 = 𝑣 𝐸𝑦 , 𝑑 𝑣 𝐸𝑦 𝑑𝑡 =− 𝐺 𝑀 𝑆 𝑦 𝐸 𝑟 𝐸 3 − 𝐺 𝑀 𝐽 𝑦 𝐸 − 𝑦 𝐽 𝑟 𝐸𝐽 3 𝑟 𝐸 = 𝑥 𝐸 2 + 𝑦 𝐸 2 , 𝑟 𝐸𝐽 = 𝑥 𝐸 − 𝑥 𝐽 2 + 𝑦 𝐸 − 𝑦 𝐽 2

성분별 운동 방정식 𝑑 𝑥 𝐽 𝑑𝑡 = 𝑣 𝐽𝑥 , 𝑑 𝑣 𝐽𝑥 𝑑𝑡 =− 𝐺 𝑀 𝑆 𝑥 𝐽 𝑟 𝐽 3 − 𝐺 𝑀 𝐸 𝑥 𝐽 − 𝑥 𝐸 𝑟 𝐽𝐸 3 𝑑 𝑦 𝐽 𝑑𝑡 = 𝑣 𝐽𝑦 , 𝑑 𝑣 𝐽𝑦 𝑑𝑡 =− 𝐺 𝑀 𝑆 𝑦 𝐽 𝑟 𝐽 3 − 𝐺 𝑀 𝐸 𝑦 𝐽 − 𝑦 𝐸 𝑟 𝐽𝐸 3 𝑟 𝐽 = 𝑥 𝐽 2 + 𝑦 𝐽 2 , 𝑟 𝐽𝐸 = 𝑥 𝐽 − 𝑥 𝐸 2 + 𝑦 𝐽 − 𝑦 𝐸 2 = 𝑟 𝐸𝐽 Motion.y2015.Week06.ThreeBody

초기조건 근일점(Perihelion): 𝑟 𝑚𝑖𝑛 =𝑎 1−𝑒 원일점(Aphelion): 𝑟 𝑚𝑎𝑥 =𝑎 1+𝑒 𝑎=1.0,(Earth) 9.58(Jupitor) 𝑒=0.106,(Earth) 0.049(Jupitor) 𝑣 𝑚𝑎𝑥 = 𝐺 𝑀 𝑆 1+𝑒 𝑎 1−𝑒 1+ 𝑀 𝐸 𝑀 𝑆 𝑣 𝑚𝑖𝑛 = 𝐺 𝑀 𝑆 1−𝑒 𝑎 1+𝑒 1+ 𝑀 𝐸 𝑀 𝑆

Modelica 코드

Modelica 코드

결과해석 목성의 질량이 실제 질량과 같으면 두 행성의 궤도는 안정적이다. 목성의 질량이 적당히 커지면 지구가 목성의 위성처럼 행동한다. 목성의 질량이 크면 지구의 궤도의 매우 불안정하다.

태양계의 공명 : 커크우드(Kirkwood) 틈과 행성 고리 행성과 태양의 거리(Titus-Bode formula) 행성 a(AU) Titus-Bode(AU) Number 수성(Mercury) 0.39 0.40 금성(Venus) 0.72 0.70 3 지구(Earth) 1.00 6 화성(Mars) 1.52 1.60 12 ???? ... 2.80 24 목성(Jupiter) 5.20 48 토성(Saturn) 9.54 10.00 96 천왕성(Uranus) 19.19 19.60 192 해왕성(Neptune) 30.06 38.80 384 명왕성(Pluto) 39.53 77.20 768 (𝑁𝑢𝑚𝑏𝑒𝑟+4)/10

Titus-Bode 식 물리적인 식이 아니라 거리의 경향을 보여준다. 목성과 토성 사이에 행성을 예측한다.(물론 없다.) 목성과 토성 사이에 많은 소행성들(asteroids)이 있다.

커크우드 틈(Kirkwood Gap) 소행성이 없는 영역이 있다. 없는 곳은 목성과 공명인 주기를 갖는다. 3:1 7:3 2:1 5:2 from Wikipedia

커크우드 틈 목성과 공명 위치에 있는 소행성은 목성의 영향으로 궤도가 큰 영향을 받는다. 태양-소행성-목성으로 시뮬레이션할 수 있다. 이 경우 소행성 사이의 상호작용은 무시한다. 원운동 가정 초기 조건 대상 반경(AU) 속도(AU/yr) 질량(kg) 소행성1 3.000 3.628 3.0× 10 21 소행성2(2:1공명) 3.276 3.471 소행성3 3.700 3.267 목성 5.200 2.755 1.9× 10 27

시뮬레이션 결과

결과 해석 공명 위치에 있는 소행성은 궤도가 목성의 영향으로 크게 변하여 궤도를 유지할 수 없게 된다. 토성 주위의 고리도 같은 현상을 보인다.

하이퍼리온(Hyperion)의 카오스 운동 spin-orbit 공명 거의 모든 행성의 달은 한 면이 항상 행성을 향하고 있다. 달의 고르지 않은 질량 분포가 에너지를 소비한다. 따라서 회전이 늦어 진다. 달의 회전이 공전주기와 일치하면 spin-orbit 공명을 통하여 회전을 유지 시킨다. 예외가 토성의 달 하이퍼리온(Hyperion) 이다.

하이퍼리온(Hyperion) 시뮬레이션 𝑚 1 𝑥 1 , 𝑦 1 카오스 운동 𝜃 𝐹 1 𝑥 2 , 𝑦 2 𝑚 2 𝐹 2 𝜔= 𝑑𝜃 𝑑𝑡 𝑥 𝑐 , 𝑦 𝑐 = 𝑚 1 𝑥 1 + 𝑚 2 𝑥 2 𝑚 1 + 𝑚 2 , 𝑚 1 𝑦+ 𝑚 2 𝑦 2 𝑚 1 + 𝑚 2

질량 중심의 미분 방정식 𝑑 2 𝑟 𝑐 𝑑 𝑡 2 =− 𝐺 𝑀 𝑆𝑎𝑡 𝑟 𝑐 𝑟 𝑐 3 =− 4 𝜋 2 𝑟 𝑐 𝑟 𝑐 3 (Hyperion단위) 𝑑 𝑣 𝑐𝑥 𝑑𝑡 =− 4 𝜋 2 𝑥 𝑐 𝑟 𝑐 3 , 𝑑 𝑥 𝑐 𝑑𝑡 = 𝑣 𝑐𝑥 𝑑 𝑣 𝑐𝑦 𝑑𝑡 =− 4 𝜋 2 𝑦 𝑐 𝑟 𝑐 3 , 𝑑 𝑦 𝑐 𝑑𝑡 = 𝑣 𝑐𝑦

운동 방정식 힘 토크 각운 동 방정식 𝐹 1 =− 𝐺 𝑀 𝑆𝑎𝑡 𝑚 1 𝑟 1 3 𝑥 1 𝑖 + 𝑦 1 𝑗 𝐹 1 =− 𝐺 𝑀 𝑆𝑎𝑡 𝑚 1 𝑟 1 3 𝑥 1 𝑖 + 𝑦 1 𝑗 𝐹 2 =− 𝐺 𝑀 𝑆𝑎𝑡 𝑚 2 𝑟 2 3 𝑥 2 𝑖 + 𝑦 2 𝑗 토크 𝜏 1 = 𝑥 1 − 𝑥 𝑐 𝑖 + 𝑦 1 − 𝑦 𝑐 𝑗 × 𝐹 1 𝜏 2 = 𝑥 2 − 𝑥 𝑐 𝑖 + 𝑦 2 − 𝑦 𝑐 𝑗 × 𝐹 2 각운 동 방정식 𝑑 𝜔 𝑑𝑡 = 𝜏 1 + 𝜏 2 𝐼 , 𝐼= 𝑚 1 𝑟 1 2 + 𝑚 2 𝑟 2 2

각 운동 방정식 𝑑𝜔 𝑑𝑡 ≈− 3𝐺 𝑀 𝑆𝑎𝑡 𝑟 𝑐 5 𝑥 𝑐 sin 𝜃 − 𝑦 𝑐 cos 𝜃 𝑥 𝑐 cos 𝜃 + 𝑦 𝑐 sin 𝜃 𝑥 𝑐 cos 𝜃 + 𝑦 𝑐 sin 𝜃 =− 12 𝜋 2 𝑟 𝑐 5 𝑥 𝑐 sin 𝜃 − 𝑦 𝑐 cos 𝜃 𝑥 𝑐 cos 𝜃 + 𝑦 𝑐 sin 𝜃

원궤도 결과

타원궤도 결과

카오스 운동 ∆𝜃= 𝜃 1 − 𝜃 2 2

결과 해석 원 궤도인 경우에는 카오스 운동이 아니다. 타원 궤도 인 경우에 카오스 운동을 한다. 시뮬레이션은 매우 단순화 한 것으로 실제 운동과는 다르다.