J2 섭동에 의한 궤도요소의 시간 변화율 - 1

J2 섭동에 의한 궤도요소(orbital elements)의 시간 변화율은 라그랑지 행성 방정식(Lagrange planetary equation)이나 가우스 행성 방정식(Gauss planetary equation)을 이용하여 계산할 수 있다. 여기서는 가우스 행성 방정식을 이용해서 계산해 보도록 하겠다.

게시글 (https://pasus.tistory.com/346)에 있는 가우스 행성 방정식은 다음과 같았다.

 

$$\begin{array}{r}\text{(1)}& & \frac{da}{dt}=\frac{2a^2}{h}e\sin \theta a_r+\frac{2a^2}{h}(1+e\cos \theta )a_{\theta }\\ & \frac{de}{dt}=\frac{h}{\mu }\sin \theta a_r+\frac{1}{h}[(r+\frac{h^2}{\mu })\cos \theta +er]a_{\theta }\\ \\ & \frac{di}{dt}=\frac{r}{h}\cos (\omega +\theta )a_h\\ \\ & \frac{d\mathrm{\Omega }}{dt}=\frac{r\sin (\omega +\theta )}{h\sin i}{a}_h\\ \\ & \frac{d\omega }{dt}=-\frac{h}{e\mu }\cos \theta a_r+\frac{1}{eh}(r+\frac{h^2}{\mu })\sin \theta a_{\theta }-\frac{r\sin (\omega +\theta )}{h\tan i}{a}_h\\ \\ & \frac{d\theta }{dt}=\frac{h}{r^2}+\frac{h}{e\mu }\cos \theta a_r-\frac{1}{eh}(r+\frac{h^2}{\mu })\sin \theta a_{\theta }\\ \\ & \frac{dh}{dt}=r{a}_{\theta }\end{array}$$

 

여기서 $a_r,a_{\theta },a_h$ 는 LVLH(local-vertical-local-horizontal) 좌표계 $\{o\}$ 로 표현된 섭동 가속도 벡터 ${\overrightarrow{a}}_p$ 의 축 성분이다.

 

$$\begin{array}{}\text{(2)}& {\mathbf{a}}_p^o=\left[\begin{array}{c}{a}_r\\ a_{\theta }\\ a_h\end{array}\right]\end{array}$$

 

한편 게시글 (https://pasus.tistory.com/349)에 의하면 J2 섭동 가속도는 다음과 같이 ECI 죄표계 $\{i\}$ 로 주어졌다.

 

$$\begin{array}{}\text{(3)}& {\mathbf{a}}_p^i=\left[\begin{array}{c}{a}_x\\ a_y\\ a_z\end{array}\right]=\frac{3}{2}{J}_2\frac{\mu }{r^2}{\left(\frac{R_e}{r}\right)}^2\left[\begin{array}{c}\frac{x}{r}(5{\left(\frac{z}{r}\right)}^2-1)\\ \frac{y}{r}(5{\left(\frac{z}{r}\right)}^2-1)\\ \frac{z}{r}(5{\left(\frac{z}{r}\right)}^2-3)\end{array}\right]\end{array}$$

 

여기서 $x,y,z$ 는 위치벡터 $\overrightarrow{r}$ 을 ECI 좌표계로 표현한 축 성분이다.

식 (1)과 (3)에 의하면 가우스 행성 방정식을 사용하려면 J2 섭동 가속도를 LVLH 좌표계로 변환해야 한다. ECI 좌표계 $\{i\}$ 에서 LVLH 좌표계 $\{o\}$ 로의 DCM $C_o^i$ 는 다음과 같다.

 

$$\begin{array}{r}\text{(4)}& C_o^i& =C(z,\mathrm{\Omega })C(x,i)C(z,\omega +\theta )\\ & =\left[\begin{array}{ccc}c\mathrm{\Omega }c(\omega +\theta )-s\mathrm{\Omega }cis(\omega +\theta )& -c\mathrm{\Omega }s(\omega +\theta )-s\mathrm{\Omega }cic(\omega +\theta )& s\mathrm{\Omega }si\\ s\mathrm{\Omega }c(\omega +\theta )+c\mathrm{\Omega }cis(\omega +\theta )& -s\mathrm{\Omega }s(\omega +\theta )+c\mathrm{\Omega }cic(\omega +\theta )& -c\mathrm{\Omega }si\\ sis(\omega +\theta )& sic(\omega +\theta )& ci\end{array}\right]\end{array}$$

 

여기서 $c=\cos ,s=\sin$ 을 의미한다.

 

 

그러면 섭동 가속도는 다음과 같이 변환할 수 있다.

 

$$\begin{array}{r}\text{(5)}& & {\mathbf{a}}_p^o={\left(C_o^i\right)}^T{\mathbf{a}}_p^i\\ & \to \left[\begin{array}{c}{a}_r\\ a_{\theta }\\ a_h\end{array}\right]=\left[\begin{array}{c}c\alpha c\varphi a_x+s\alpha c\varphi a_y+s\varphi a_z\\ -s\alpha a_x+c\alpha a_y\\ -c\alpha s\varphi a_x-s\alpha s\varphi a_y+c\varphi a_z\end{array}\right]\end{array}$$

 

하지만 여기서 섭동 가속도를 좌표변환하기 전에 먼저 식 (3)의 섭동 가속도식에 있는 ECI 좌표계와 관련된 변수 $x,y,z$ 를 LVLH 좌표계 변수와 관련된 궤도요소의 함수로 바꿔야 한다. 이를 위해서 임시 좌표계 $\{t\}$ 를 도입한다. 임시 좌표계는 다음과 같이 ECI 좌표계의 z축을 중심으로 $\alpha$ 만큼 회전한 후 다시 y축으로 $-\varphi$ 만큼 회전하여 얻어진 좌표계다. 따라서 임시 좌표계의 x축은 LVLH 좌표계의 x축과 일치한다.

 

 

그러면 ECI 좌표계에서 임시 좌표계로의 DCM $C_t^i$ 는 다음과 같다.

 

$$\begin{array}{r}\text{(6)}& C_t^i& =C(z,\alpha )C(y,-\varphi )\\ & =\left[\begin{array}{ccc}c\alpha c\varphi & -s\alpha & -c\alpha s\varphi \\ s\alpha c\varphi & c\alpha & -s\alpha s\varphi \\ s\varphi & 0& c\varphi \end{array}\right]\end{array}$$

 

좌표변환을 이용하면 위치벡터 $\overrightarrow{r}$ 은 다음과 같이 각각 임시 좌표계, LVLH 좌표계 및 ECI 좌표계로 표현할 수 있다.

 

$$\begin{array}{}\text{(7)}& {\mathbf{r}}^o={\mathbf{r}}^t=\left[\begin{array}{c}r\\ 0\\ 0\end{array}\right],{\mathbf{r}}^i=C_o^i{\mathbf{r}}^o=C_t^i{\mathbf{r}}^t\end{array}$$

 

따라서 식 (4), (6), (7)을 이용하면 다음과 같은 관계식을 얻을 수 있다.

 

$$\begin{array}{r}\text{(8)}& \frac{x}{r}& =\cos \alpha \cos \varphi \\ & =\cos \mathrm{\Omega }\cos (\omega +\theta )-\sin \mathrm{\Omega }\cos i\sin (\omega +\theta )\\ \\ \frac{y}{r}& =\sin \alpha \cos \varphi \\ \\ & =\sin \mathrm{\Omega }\cos (\omega +\theta )+\cos \mathrm{\Omega }\cos i\sin (\omega +\theta )\\ \\ \frac{z}{r}& =\sin \varphi =\sin i\sin (\omega +\theta )\end{array}$$

 

식 (8)을 식 (3)에 대입하면 J2 섭동 가속도는 다음과 같이 된다.

 

$$\begin{array}{}\text{(9)}& {\mathbf{a}}_p^i=\frac{3}{2}{J}_2\frac{\mu }{r^2}{\left(\frac{R_e}{r}\right)}^2\left[\begin{array}{c}\cos \alpha \cos \varphi (5{\sin }^2\varphi -1)\\ \sin \alpha \cos \varphi (5{\sin }^2\varphi -1)\\ \sin \varphi (5{\sin }^2\varphi -3)\end{array}\right]\end{array}$$

 

 

 

이제 식 (9)를 (5)에 대입하면 다음과 같이 LVLH 좌표계로 표현된 J2 섭동 가속도의 축 성분을 구할 수 있다.

 

$$\begin{array}{rl}\text{(10)}& {\tilde{a}}_r& =(c\mathrm{\Omega }c(\omega +\theta )-s\mathrm{\Omega }cis(\omega +\theta ))a_x& +(s\mathrm{\Omega }c(\omega +\theta )+c\mathrm{\Omega }cis(\omega +\theta ))a_y\\ & +(sis(\omega +\theta ))a_z\\ \\ & =\cos \alpha \cos \varphi a_x+\sin \alpha \cos \varphi a_y+\sin \varphi a_z\\ \\ & ={\cos }^2\alpha {\cos }^2\varphi (5{\sin }^2\varphi -1)+{\sin }^2\alpha {\cos }^2\varphi (5{\sin }^2\varphi -1)\\ & +{\sin }^2\varphi (5{\sin }^2\varphi -3)\\ \\ & =(5{\sin }^2\varphi -1)-2{\sin }^2\varphi \\ \\ & =3{\sin }^{2}i{\sin }^2(\omega +\theta )-1\\ \\ \\ \text{(11)}& {\tilde{a}}_{\theta }& =(-c\mathrm{\Omega }s(\omega +\theta )-s\mathrm{\Omega }cic(\omega +\theta ))a_x& +(-s\mathrm{\Omega }s(\omega +\theta )+c\mathrm{\Omega }cic(\omega +\theta ))a_y\\ & +(sic(\omega +\theta ))a_z\\ \\ & =\left[\begin{array}{c}-c\mathrm{\Omega }s(\omega +\theta )-s\mathrm{\Omega }cic(\omega +\theta ))\cos \alpha \cos \varphi \\ +(-s\mathrm{\Omega }s(\omega +\theta )+c\mathrm{\Omega }cic(\omega +\theta ))\sin \alpha \cos \varphi \\ +sic(\omega +\theta )\sin \varphi \end{array}\right](5{\sin }^2\varphi -1)\\ & -2\sin i\cos (\omega +\theta )\sin \varphi \\ \\ & =-2{\sin }^{2}i\sin (\omega +\theta )\cos (\omega +\theta )\\ \\ & =-{\sin }^{2}i\sin 2(\omega +\theta )\\ \\ \\ \text{(12)}& {\tilde{a}}_h& =s\mathrm{\Omega }si{a}_x-c\mathrm{\Omega }si{a}_y+ci{a}_z\\ & =\sin \mathrm{\Omega }\sin i\cos \alpha \cos \varphi (5{\sin }^2\varphi -1)\\ & -\cos \mathrm{\Omega }\sin i\sin \alpha \cos \varphi (5{\sin }^2\varphi -1)+\cos i\sin \varphi (5{\sin }^2\varphi -3)\\ \\ & =\left[\begin{array}{c}\sin \mathrm{\Omega }\sin i\cos \alpha \cos \varphi \\ -\cos \mathrm{\Omega }\sin i\sin \alpha \cos \varphi \\ +\cos i\sin \varphi \end{array}\right](5{\sin }^2\varphi -1)-2\cos i\sin \varphi \\ \\ & =-2\cos i\sin \varphi \\ \\ & =-\sin 2i\sin (\omega +\theta )\end{array}$$

 

여기서 ${\tilde{a}}_r,{\tilde{a}}_{\theta },{\tilde{a}}_h$ 에 $\frac{3}{2}{J}_2\frac{\mu }{r^2}{\left(\frac{R_e}{r}\right)}^2$ 을 곱하면 $a_r,a_{\theta },a_h$ 을 얻을 수 있다. 또한 식 (8)을 이용하면 식 (11)과 (12)에 있는 대괄호 항은 $0$ 이 됨을 보일 수 있다.

촤종적으로 가우스 행성 방정식 (1)에 식 (11), (12), (13)을 대입하고 전개하면 J2 섭동에 의한 궤도요소의 시간 변화율은 다음과 같다.

 

$$\begin{array}{r}\text{(13)}& \frac{da}{dt}& =3J_2\frac{a^2\mu R_e^2}{h{r}^4}\left[\begin{array}{c}e\sin \theta (3{\sin }^{2}i{\sin }^2(\omega +\theta )-1)\\ -(1+e\cos \theta ){\sin }^{2}i\sin 2(\omega +\theta )\end{array}\right]\\ \frac{de}{dt}& =\frac{3}{2}{J}_2\frac{\mu R_e^2}{h{r}^3}\left[\begin{array}{c}\frac{h^2}{\mu r}\sin \theta (3{\sin }^{2}i{\sin }^2(\omega +\theta )-1)\\ -[(2+e\cos \theta )\cos \theta +e]{\sin }^{2}i\sin 2(\omega +\theta )\end{array}\right]\\ \\ \frac{di}{dt}& =\frac{3}{2}{J}_2\frac{\mu R_e^2}{eh{r}^3}\left[\begin{array}{c}\frac{h^2}{\mu r}\cos \theta (1-3{\sin }^{2}i{\sin }^2(\omega +\theta ))\\ -(2+e\cos \theta )\sin \theta {\sin }^{2}i\sin 2(\omega +\theta )\\ +2e{\cos }^{2}i{\sin }^2(\omega +\theta )\end{array}\right]\\ \\ \frac{d\mathrm{\Omega }}{dt}& =-3J_2\frac{\mu R_e^2}{h{r}^3}\cos i{\sin }^2(\omega +\theta )\\ \\ \frac{d\omega }{dt}& =\frac{3}{2}{J}_2\frac{\mu R_e^2}{eh{r}^3}\left[\begin{array}{c}\frac{h^2}{\mu r}\cos \theta (1-3{\sin }^{2}i{\sin }^2(\omega +\theta ))\\ -(2+e\cos \theta )\sin \theta {\sin }^{2}i\sin 2(\omega +\theta )\\ +2e{\cos }^{2}i{\sin }^2(\omega +\theta )\end{array}\right]\\ \\ \frac{d\theta }{dt}& =\frac{h}{r^2}+\frac{3}{2}{J}_2\frac{\mu R_e^2}{eh{r}^3}\left[\begin{array}{c}\frac{h^2}{\mu r}\cos \theta (3{\sin }^{2}i{\sin }^2(\omega +\theta )-1)\\ +(2+e\cos \theta )\sin \theta {\sin }^{2}i\sin 2(\omega +\theta )\end{array}\right]\\ \\ \frac{dh}{dt}& =-\frac{3}{2}{J}_2\frac{\mu R_e^2}{r^3}{\sin }^{2}i\sin 2(\omega +\theta )\end{array}$$