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

지구의 운동 1차 방정식으로 변환 Modelica 기술 𝑑 𝑣 𝑥 𝑑𝑡 =− 4 𝜋 2 𝑥 𝑟 3 , 𝑑𝑥 𝑑𝑡 = 𝑣 𝑥 𝑑 𝑣 𝑥 𝑑𝑡 =− 4 𝜋 2 𝑥 𝑟 3 , 𝑑𝑥 𝑑𝑡 = 𝑣 𝑥 𝑑 𝑣 𝑦 𝑑𝑡 =− 4 𝜋 2 𝑦 𝑟 3 , 𝑑𝑦 𝑑𝑡 = 𝑣 𝑦 Modelica 기술 𝑑𝑒𝑟 𝑣𝑥 =−4∗𝑃 𝐼 2 ∗𝑥/𝑟^3; 𝑑𝑒𝑟 𝑥 =𝑣𝑥; 𝑑𝑒𝑟 𝑣𝑦 =−4∗𝑃 𝐼 2 ∗𝑦/𝑟^3; 𝑑𝑒𝑟 𝑦 =𝑣𝑦;

불필요한 Package/Class 닫기 Motion 패키지가 사라졌음

Package 생성 패키지 명 Motion 패키지 선택 새로운 패키지 Motion이 root(global) 에 생성됨

Package 생성 패키지 명 y2016 패키지 선택 상위 패키지 새로운 패키지 y2016가 Motion 아래에 생성됨

Package 생성 Week06 Week06

Package 저장 D:\lec_hwlee\motion\y2016\week06 Motion.y2016.Week06.mo

클래스 생성 새로운 클래스 Earth가 Motion.y2016.Week06 아래에 생성됨 Earth Earth

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

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

시뮬레이션 조건 설정 Earth

시뮬레이션 조건 설정 시뮬레이션 시작시간 종료시간 설정 미분방정식 Solver 선택 오차 한계 y2016.Week06.Earth 시뮬레이션 시작시간 종료시간 설정 미분방정식 Solver 선택 dassl 오차 한계

시뮬레이션 조건 설정 결과 출력 형식 mat : 이진 형식 출력, MATLAB, Octave 에서 사용 가능 y2016.Week06.Earth 결과 출력 형식 mat : 이진 형식 출력, MATLAB, Octave 에서 사용 가능 plt : 일반 텍스트 출력 csv : 자료를 콤마로 구분하여 저장 empty : 출력하지 않음 적분 구간의 개수 ∆𝑡= 종료시간 −시작시간 인터벌 수 ∆𝑡= 10−0 10000 =0.001초

시뮬레이션 실행 콤파일하고 시뮬레이션한 로그를 보여준다.

시뮬레이션 실행

결과 그래프 꾸미기 Setup을 통하여 그래프의 형식을 조정할 수 있다.

결과 그래프 꾸미기 그래프 타이틀 y축 타이틀 x축 타이틀 범례 선색 선굵기

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

결과 보이기

다른 파라메터 값의 결과 보기 파라메터 값을 수정하고 re-simulate을 한다. 오른버튼 클릭

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

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

다른 Solver 적용 y2016.Week06.Earth y2016.Week06.Earth Runge-Kutta 선택

RK 결과

역 제곱 법칙 역 제곱 법칙은 정확한 것인가? 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 일반상대론: 46 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 코드

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

삼체 문제와 목성의 지구에 대한 효과 𝐹 𝐸𝐽 𝐹 𝐸𝑆 𝐹 𝐸𝐽 지구 목성 𝐹 𝐸𝑆 태양 𝐹 𝐸𝐽 =− 𝐺 𝑀 𝐸 𝑀 𝐽 𝑟 𝐸𝐽 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 코드 Motion.y2016.Week06.ThreeBody

Modelica 코드

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

커크우드 틈(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) 시뮬레이션 𝑚 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 𝜃

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

도전해 보기 초기조건이 다른 행성의 운동을 한 그림으로 그리기 연습문제 4.1~4.7 풀이하기 연습문제 4.8, 4.9 풀이하기 연습문제 4.10, 4.11 풀이하기 삼체 문제에서 목성이 있는 경우와 없는 경우에 지구의 궤도의 차이를 계산하는 클래스를 만드시오.

도전해 보기 연습문제 4.12~4.16 풀이하기 연습문제 4.17, 4.18 풀이하기 연습문제 4.19, 4.20 풀이하기