가우스 변분 방정식 (Gauss Variational Equation)

라그랑지 행성 방정식은 섭동력이 보존력(conservative force)이어야 한다는 조건이 있었다. 하지만 섭동력이 보존력이 아닌 경우도 많다. 예를 들면 대기 항력, 제어 추력, 태양 복사 압력 등이다. 특히 섭동력이 제어 추력인 경우, 이 힘이 궤도요소에 어떤 영향을 미치는지를 직접적으로 이해하는 것은 제어기 설계에 있어서 매우 중요하다.

가우스 변분 방정식(Gauss variational equation)은 임의의 섭동력으로 인한 궤도요소의 시간 변화율을 힘의 관점에서 명시적으로 표현하기 때문에 섭동력이 비보존력인 경우에 특히 유용하다. 더구나 보존력인 경우에도 힘을 포텐셜 함수의 그래디언트로 표현할 수 있기 때문에 적용 가능하다. 라그랑지 행성 방정식을 유도할 때는 특별한 좌표계를 언급하지 않은 채 관성좌표계를 사용했지만 여기서는 LVLH (local-vertical-local-horizontal) 좌표계를 이용한다. LVLH 좌표계는 가우스 변분 방정싱을 가장 쉽게 유도할 수 있는 좌표계이다. 좌표계는 $\{o\}$ 로 표시하며 우주비행체의 중심에 원점이 있고 x( ${\hat{o}}_1$)축은 지구 중심으로부터 비행체의 위치 벡터 방향, z(${\hat{o}}_3$) 축은 각운동량 벡터 방향, y(${\hat{o}}_2$) 축은 오른손 법칙에 의해서 정해진다 (https://pasus.tistory.com/238).

 

 

임의의 섭동 가속도를 포함하는 이체문제 운동방정식은 다음과 같다.

 

$$\begin{array}{}\text{(1)}& \frac{{}^i{d}^2\overrightarrow{r}}{d{t}^2}+\frac{\mu }{r^3}\overrightarrow{r}={\overrightarrow{a}}_p\end{array}$$

 

여기서 ${\overrightarrow{a}}_p$ 는 섭동 가속도이다. 이체문제에서는 각운동량 벡터 $\overrightarrow{h}$, 이심율 벡터(eccentricity vector) $\overrightarrow{e}$, 승교선 벡터(ascending node vector) $\overrightarrow{n}$ 는 모두 관성좌표계에 대해서 상수벡터이지만 (https://pasus.tistory.com/287), 섭동 가속도가 있는 경우에는 그렇지 않다.

먼저 임의의 섭동력에 의해서 각운동량 벡터가 어떻게 변화하는지 알아보자. 각운동량 벡터는 $\overrightarrow{h}=\overrightarrow{r}\times \overrightarrow{v}$ 이므로 식 (1)을 이용하여 시간 미분하면 다음과 같다.

 

$$\begin{array}{r}\text{(2)}& \frac{{}^{i}d\overrightarrow{h}}{dt}& =\frac{{}^{i}d\overrightarrow{r}}{dt}\times \overrightarrow{v}+\overrightarrow{r}\times \frac{{}^{i}d\overrightarrow{v}}{dt}\\ & =\overrightarrow{v}\times \overrightarrow{v}+\overrightarrow{r}\times (-\frac{\mu }{r^3}\overrightarrow{r}+{\overrightarrow{a}}_p)\\ \\ & =\overrightarrow{r}\times {\overrightarrow{a}}_p\end{array}$$

 

이심율 벡터는 $\overrightarrow{e}=\frac{1}{\mu }(\overrightarrow{v}\times \overrightarrow{h}-\mu \frac{\overrightarrow{r}}{r})$ 이므로, 미분하면 다음과 같다.

 

$$\begin{array}{r}\text{(3)}& \frac{{}^{i}d\overrightarrow{e}}{dt}& =\frac{1}{\mu }(\frac{{}^{i}d\overrightarrow{v}}{dt}\times \overrightarrow{h}+\overrightarrow{v}\times \frac{{}^{i}d\overrightarrow{h}}{dt})-\frac{{}^{i}d}{dt}\left(\frac{\overrightarrow{r}}{r}\right)\\ & =\frac{1}{\mu }((-\frac{\mu }{r^3}\overrightarrow{r}+{\overrightarrow{a}}_p)\times \overrightarrow{h}+\overrightarrow{v}\times (\overrightarrow{r}\times {\overrightarrow{a}}_p))-\frac{{}^{i}d}{dt}\left(\frac{\overrightarrow{r}}{r}\right)\\ \\ & =\frac{1}{\mu }({\overrightarrow{a}}_p\times \overrightarrow{h}+\overrightarrow{v}\times (\overrightarrow{r}\times {\overrightarrow{a}}_p))-[\frac{{}^{i}d}{dt}\left(\frac{\overrightarrow{r}}{r}\right)+\frac{1}{r^3}\overrightarrow{r}\times \overrightarrow{h}]\end{array}$$

 

그런데 여기서

 

$$\begin{array}{rl}\frac{{}^{i}d}{dt}\left(\frac{\overrightarrow{r}}{r}\right)+\frac{1}{r^3}\overrightarrow{r}\times \overrightarrow{h}& =\frac{\overrightarrow{v}}{r}-\frac{\overrightarrow{r}}{r^3}(\overrightarrow{r}\cdot \overrightarrow{v})+\frac{1}{r^3}\overrightarrow{r}\times (\overrightarrow{r}\times \overrightarrow{v})\\ \\ & =\frac{\overrightarrow{v}}{r}-\frac{\overrightarrow{r}}{r^3}(\overrightarrow{r}\cdot \overrightarrow{v})+\frac{1}{r^3}[\overrightarrow{r}(\overrightarrow{r}\cdot \overrightarrow{v})-\overrightarrow{v}(\overrightarrow{r}\cdot \overrightarrow{r})]\\ \\ & =\frac{\overrightarrow{v}}{r}+\frac{1}{r^3}[-r^2\overrightarrow{v}]\\ \\ & =0\end{array}$$

 

이므로 식 (3)은 다음과 같이 된다.

 

$$\begin{array}{}\text{(4)}& \frac{{}^{i}d\overrightarrow{e}}{dt}=\frac{1}{\mu }({\overrightarrow{a}}_p\times \overrightarrow{h}+\overrightarrow{v}\times (\overrightarrow{r}\times {\overrightarrow{a}}_p))\end{array}$$

 

승교선 벡터는 $\overrightarrow{n}={\hat{i}}_3\times \overrightarrow{h}$ 이므로 미분하면 다음과 같다.

 

$$\begin{array}{}\text{(5)}& \frac{{}^{i}d\overrightarrow{n}}{dt}={\hat{i}}_3\times \frac{{}^{i}d\overrightarrow{h}}{dt}={\hat{i}}_3\times (\overrightarrow{r}\times {\overrightarrow{a}}_p)\end{array}$$

 

BKE에 의하면 각운동량 벡터, 이심율 벡터, 승교선 벡터를 관성좌표계에서의 미분과 LVLH 좌표계에서의 미분 간에는 다음 관계식이 성립한다 (https://pasus.tistory.com/121).

 

$$\begin{array}{r}\text{(6)}& & \frac{{}^{i}d\overrightarrow{h}}{dt}=\frac{{}^{o}d\overrightarrow{h}}{dt}{+}^i{\overrightarrow{\omega }}^o\times \overrightarrow{h}\\ & \frac{{}^{i}d\overrightarrow{e}}{dt}=\frac{{}^{o}d\overrightarrow{e}}{dt}{+}^i{\overrightarrow{\omega }}^o\times \overrightarrow{e}\\ \\ & \frac{{}^{i}d\overrightarrow{n}}{dt}=\frac{{}^{o}d\overrightarrow{n}}{dt}{+}^i{\overrightarrow{\omega }}^o\times \overrightarrow{n}\end{array}$$

 

여기서 ${}^i{\overrightarrow{\omega }}^o$ 는 관성좌표계를 기준으로 한 LVLH 좌표계의 각속도 벡터로서 다음과 같다 (https://pasus.tistory.com/120).

 

$$\begin{array}{}\text{(7)}& {}^i{\overrightarrow{\omega }}^o=\dot{\mathrm{\Omega }}{\hat{i}}_3+\dot{i}\hat{n}+(\dot{\omega }+\dot{\theta }){\hat{o}}_3\end{array}$$

 

여기서 $\hat{n}=\frac{\overrightarrow{n}}{n}$ 이다.

 

 

한편 LVLH 좌표계 정의에 의하면 $\overrightarrow{r}=r{\hat{o}}_1$, $\overrightarrow{h}=h{\hat{o}}_3$ 이며, 섭동 가속도도 LVLH 좌표계로 다음과 같이 표현할 수 있다.

 

$$\begin{array}{}\text{(8)}& {\overrightarrow{a}}_p=a_r{\hat{o}}_1+a_{\theta }{\hat{o}}_2+a_h{\hat{o}}_3\end{array}$$

 

 

 

이제 모든 벡터를 LVLH 좌표계로 표현해 보자.

 

$$\begin{array}{r}\text{(9)}& {\mathbf{r}}^o& =\left[\begin{array}{c}r\\ 0\\ 0\end{array}\right],{\mathbf{h}}^o=\left[\begin{array}{c}0\\ 0h\end{array}\right],{\mathbf{a}}_p^o=\left[\begin{array}{c}{a}_r\\ a_{\theta }\\ a_h\end{array}\right]\\ {\omega }_{io}^o& =C_i^o\left[\begin{array}{c}0\\ 0\\ \dot{\mathrm{\Omega }}\end{array}\right]+C_n^o\left[\begin{array}{c}\dot{i}\\ 0\\ 0\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ \dot{\omega }+\dot{\theta }\end{array}\right]\\ \\ & =(C(z,\mathrm{\Omega })C(x,i)C(z,\omega +\theta ){)}^T\left[\begin{array}{c}0\\ 0\\ \dot{\mathrm{\Omega }}\end{array}\right]\\ \\ & +(C(x,i)C(z,\omega +\theta ){)}^T\left[\begin{array}{c}\dot{i}\\ 0\\ 0\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ \dot{\omega }+\dot{\theta }\end{array}\right]\\ \\ & =\left[\begin{array}{c}\dot{\mathrm{\Omega }}\sin (\omega +\theta )\sin i+\dot{i}\cos (\omega +theta)\\ \dot{\mathrm{\Omega }}\cos (\omega +\theta )\sin i-\dot{i}\sin (\omega +\theta )\\ \dot{\mathrm{\Omega }}\cos i+\dot{\omega }+\dot{\theta }\end{array}\right]\\ \\ & {\mathbf{n}}^o=\left[\begin{array}{c}h\cos (\omega +\theta )\sin i\\ -h\sin (\omega +\theta )\sin i\\ 0\end{array}\right]\end{array}$$

 

여기서 근점편각 $\omega$ 와 실제비행각 $\theta$ 의 합인 $\omega +\theta$ 를 argument of latitude 라고 한다.

속도벡터의 경우 라그랑지 제약조건(Lagrange’s constraint) 또는 오스큘레이션 제약조건(osculation constraint)에 의해서 교란되지 않은 운동의 속도벡터와 동일해야 하므로,

 

$$\begin{array}{}\text{(10)}& \overrightarrow{v}=\frac{{}^{i}d\overrightarrow{r}}{dt}=\frac{{}^{o}d\overrightarrow{r}}{dt}{+}^i{\overrightarrow{\omega }}_{osc}^o\times \overrightarrow{r}\end{array}$$

 

이어야 한다. 여기서 ${}^i{\overrightarrow{\omega }}_{osc}^o=\dot{\theta }{\hat{o}}_3$ 이다. 따라서

 

$$\begin{array}{}\text{(11)}& {\mathbf{v}}^o=\left[\begin{array}{c}\dot{r}\\ r\dot{\theta }\\ 0\end{array}\right]=\left[\begin{array}{c}{v}_r\\ v_{\theta }\\ 0\end{array}\right]\end{array}$$

 

이다. 식 (11)과 다음 관계식을 이용하면,

 

$$\begin{array}{}\text{(12)}& r=\frac{\frac{h^2}{\mu }}{1+e\cos \theta },h=r^2\dot{\theta },\dot{r}=\frac{\mu }{h}e\sin \theta \end{array}$$

 

이심율 벡터는 다음과 같이 표현할 수 있다.

 

$$\begin{array}{}\text{(13)}& {\mathbf{e}}^o=\frac{1}{\mu }\left[\begin{array}{c}h{v}_{\theta }-\mu \\ -h{v}_r\\ 0\end{array}\right]=\left[\begin{array}{c}e\cos \theta \\ -e\sin \theta \\ 0\end{array}\right]\end{array}$$

 

이제 궤도요소의 시간 변화율을 계산할 차례다. 식 (2)와 (6)에 의하면,

 

$$\begin{array}{}\text{(14)}& \frac{{}^{i}d\overrightarrow{h}}{dt}=\frac{{}^{o}d\overrightarrow{h}}{dt}{+}^i{\overrightarrow{\omega }}^o\times \overrightarrow{h}=\overrightarrow{r}\times {\overrightarrow{a}}_p\end{array}$$

 

이다. 위 식을 식 (9)를 이용하여 LVLH 좌표계로 표현하면 다음과 같다.

 

$$\begin{array}{}\text{(15)}& \left[\begin{array}{c}0\\ 0\\ \dot{h}\end{array}\right]+\left[\begin{array}{c}h(\dot{\mathrm{\Omega }}\cos (\omega +\theta )\sin i-\dot{i}\sin (\omega +\theta ))\\ -h(\dot{\mathrm{\Omega }}\sin (\omega +\theta )\sin i+\dot{i}\cos (\omega +\theta ))\\ 0\end{array}\right]=\left[\begin{array}{c}0\\ -r{a}_h\\ r{a}_{\theta }\end{array}\right]\end{array}$$

 

위 식을 이용하면 다음과 같은 식을 얻을 수 있다.

 

$$\begin{array}{r}\text{(16)}& & \frac{d\mathrm{\Omega }}{dt}=\frac{r\sin (\omega +\theta )}{h\sin i}{a}_h\\ & \frac{di}{dt}=\frac{r}{h}\cos (\omega +\theta )a_h\\ \\ & \frac{dh}{dt}=r{a}_{\theta }\end{array}$$

 

또한 식 (4)와 (6)에 의하면,

 

$$\begin{array}{}\text{(17)}& \frac{{}^{i}d\overrightarrow{e}}{dt}=\frac{{}^{o}d\overrightarrow{e}}{dt}{+}^i{\overrightarrow{\omega }}^o\times \overrightarrow{e}=\frac{1}{\mu }({\overrightarrow{a}}_p\times \overrightarrow{h}+\overrightarrow{v}\times (\overrightarrow{r}\times {\overrightarrow{a}}_p))\end{array}$$

 

이다. 위 식을 식 (9)와 (13)을 이용하여 LVLH 좌표계로 표현하면 다음과 같다.

 

$$\begin{array}{r}\text{(18)}& & \left[\begin{array}{c}\dot{e}\cos \theta -e\dot{\theta }\sin \theta \\ -\dot{e}\sin \theta -e\dot{\theta }\cos \theta \\ 0\end{array}\right]+\left[\begin{array}{c}e\sin \theta (\dot{\mathrm{\Omega }}\cos i+\dot{\omega }+\dot{\theta })\\ e\cos \theta (\dot{\mathrm{\Omega }}\cos i+\dot{\omega }+\dot{\theta })\\ e(-\dot{\mathrm{\Omega }}̇\sin i\cos \omega +\dot{i}\sin \omega )\end{array}\right]\\ & =\frac{1}{\mu }(\left[\begin{array}{c}{a}_{\theta }h\\ -a_{r}h\\ 0\end{array}\right]+\left[\begin{array}{c}r{v}_{\theta }{a}_{\theta }\\ -r{v}_r{a}_{\theta }\\ -r{v}_r{a}_h\end{array}\right])\end{array}$$

 

위 식과 식 (16)을 이용하면 다음과 같은 수식을 얻을 수 있다.

 

$$\begin{array}{r}\text{(19)}& & \frac{de}{dt}=\frac{h}{\mu }\sin \theta a_r+\frac{1}{h}[(r+\frac{h^2}{\mu })\cos \theta +er]a_{\theta }\\ & \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 }\end{array}$$

 

장반경의 변화율은 장반경 식 $a=\frac{h^2}{\mu (1-e^2)}$ 을 미분하고 $\dot{e},\dot{h}$ 을 대입하여 구할 수 있다.

 

$$\begin{array}{r}\text{(20)}& \frac{da}{dt}& =\frac{2h\dot{h}}{\mu (1-e^2)}+\frac{h^{2}2e\dot{e}}{\mu (1-e^2{)}^2}\\ & =\frac{2a^2}{h}e\sin \theta a_r+\frac{2a^2}{h}(1+e\cos \theta )a_{\theta }\end{array}$$

 

 

 

지금까지 유도한 궤도요소의 변화율을 정리하면 다음과 같다.

 

$$\begin{array}{r}\text{(21)}& & \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}$$

 

식 (21)을 가우스 변분 방정식 또는 가우스 행성 방정식(Gauss planetary equation)이라고 한다.

라그랑지 행성 방정식과 마찬가지로 위 식은 경사각이 $i=0^0$ 일 때, 이심율이 $e=0$ 일 때 특이값을 갖는다.

식 (21)을 이용하면 시간 $t_0$ 에서 6개의 궤도요소 $(a_0,e_0,i_0,{\mathrm{\Omega }}_0,{\omega }_0,{\theta }_0)$ 와 섭동 가속도 ${\overrightarrow{a}}_p$ 의 함수 형태가 주어질 때 6개의 행성 방정식을 수치적으로 적분하여 이후 시간의 궤도요소를 구할 수 있다.

또한 식 (21)은 제어 관점에서 제어력이 궤도요소에 어떤 영향을 미치는지를 직접적으로 이해하는 데 도움이 된다. 예를 들어, RAAN $\mathrm{\Omega }$ 와 경사각 $i$ 의 미분 방정식을 살펴보면 RAAN 보정을 하는 데 가장 효율적인 지점은 $\sin (\omega +\theta )$ 가 최대가 되는 곳이고, 궤도 경사를 조정하는 데 가장 효율적인 지점은 $\cos (\omega +\theta )$ 가 최대가 되는 적도를 통과하는 순간이라는 것을 알 수 있다.