[CR3BP] 주기궤도 (Periodic Orbit)의 조건

라그랑지 포인트 L1, L2 및 L3에서 선형화 운동방정식의 해석 결과, 초기값을 잘 설정한다면 주기궤도(periodic orbit)가 형성될 수 있다는 것을 알았다 (https://pasus.tistory.com/273). 하지만 선형화 운동방정식은 라그랑지 포인트에서 가까운 영역에서만 유효하기 때문에 보다 넓은 범위에서도 주기궤도를 만들 수 있는지는 더 분석해 봐야 한다.

다시 CR3BP의 무차원화된 비선형 운동방정식으로 돌아가 보자. (https://pasus.tistory.com/147).

 

\[ \begin{align} & \ddot{x}-2 \dot{y}-x= - \frac{(1-\mu)(x+\mu)}{r_1^3 }- \frac{\mu (x+\mu-1)}{r_2^3} \tag{1} \\ \\ & \ddot{y}+2 \dot{x}-y= - \frac{(1-\mu)y}{r_1^3 }- \frac{\mu y}{ r_2^3 } \\ \\ & \ddot{z}=- \frac{(1-\mu)z}{r_1^3 }- \frac{\mu z}{ r_2^3 } \end{align} \]

여기서

\[ \begin{align} & r_1= \sqrt{(x+\mu)^2+y^2+z^2 } \tag{2} \\ \\ & r_2= \sqrt{(x+\mu-1)^2+y^2+z^2 } \end{align} \]

 

이다.

식 (1)은 \(y \to -y, \ \ t \to -t\) 의 변환에 대해서 불변(invariance)이다. 즉 해당 변환에 대해서 미분 방정식이 동일하다는 얘기다. 확인해 보기 위해서 \(y=-Y, \ t=-\tau\) 로 놓자. 그러면

 

\[ \begin{align} & \dot{x}= \frac{dx}{dt}= \frac{dx}{-d \tau}= - x', \ \ \ \ \dot{y}=\frac{dy}{dt}= \frac{-dY}{-d\tau}=Y' \tag{3} \\ \\ & \ddot{x}= \frac{d^2 x}{dt^2 }= \frac{d^2x}{d\tau^2 }=x'' \\ \\ & \ddot{y}= \frac{d}{dt} \left( \frac{dy}{dt} \right) = \frac{d}{-d \tau} \left( \frac{dY}{d\tau} \right)= - \frac{d^2 Y}{d \tau^2 }= -Y'' \\ \\ & \ddot{z}= \frac{d^2 z}{dt^2 } = \frac{d^2z}{d \tau^2 }=z'' \end{align} \]

 

이 된다. 식 (3)을 식 (2)에 대입하면 동일한 \(r_1\) 과 \(r_2\) 값이 얻어지고 식 (1)에 대입하면 미분방정식은 다음과 같이 된다.

 

\[ \begin{align} & x''-2 Y'-x= - \frac{(1-\mu)(x+\mu)}{r_1^3 }- \frac{\mu (x+\mu-1)}{r_2^3} \tag{4} \\ \\ - & Y''-2 x' +Y= \frac{(1-\mu) Y}{r_1^3 }+ \frac{\mu Y}{ r_2^3 } \\ \\ & z''=- \frac{(1-\mu)z}{r_1^3 }- \frac{\mu z}{ r_2^3 } \end{align} \]

 

식 (4)와 (1)을 비교해면 동일한 식이라는 것을 알 수 있다. 즉 식 (1)은 \(y \to -y, \ \ t=-t\) 의 변환에 대해서 불변이다. 따라서 만약 식 (1)의 해(solution)가 \( (x(t), \ y(t), \ z(t), \ \dot{x}(t), \ \dot{y}(t), \ \dot{z}(t)) \) 이라면 \( (x(-t), \ -y(-t), \ z(-t), \ -\dot{x}(-t), \ \dot{y}(-t), \ -\dot{z}(-t))\) 도 해가 된다. 여기서 주목할 점은 두 해는 \((x-z)\) 평편에 대해서 대칭이라는 것이다.

 

 

따라서 만약 \((x-z)\) 평편에서 수직으로 출발하고 일정 시간이 지난 후 다시 \((x-z)\) 평편에 수직으로 도착하도록 초기조건을 설정한다면 다음 그림과 같이 \((x-z)\) 평면에 대칭인 주기궤도(periodic orbit)를 만들 수 있을 것이다.

 

 

궤적이 \((x-z)\) 평편에서 수직으로 출발해야 하므로 주기궤도를 위한 초기조건은 다음과 같아야 한다.

 

\[ \mathbf{x}(t_0 )= \begin{bmatrix} x(t_0 ) \\ y(t_0 ) \\ z(t_0 ) \\ \dot{x}(t_0 ) \\ \dot{y}(t_0 ) \\ \dot{z}(t_0 ) \end{bmatrix}= \begin{bmatrix} x_0 \\ 0 \\ z_0 \\ 0 \\ v_{y0} \\ 0 \end{bmatrix} \tag{5} \]

 

또한 주기가 \(T\) 인 주기 궤도를 위한 조건은 다음과 같아야 한다.

 

\[ \mathbf{x}(T/2 )= \begin{bmatrix} x(T/2 ) \\ 0 \\ z(T/2 ) \\ 0 \\ v_y(T/2) \\ 0 \end{bmatrix} \tag{6} \]

 

따라서 라그랑지 포인트 L1, L2 및 L3에서 주기궤도를 설계하기 위해서는 주기궤도를 위한 조건식 (6)이 성립하도록 초기조건 (5)를 설정하면 된다.

하지만 상태벡터 \(\mathbf{x}(t)\) 에 대한 초기조건의 추측값 \(\mathbf{x}(t_0 )\) 가 단번에 식 (6)의 형태대로 궤적을 생성할 가능성이 거의 없기 때문에, 궤적이 \((x-z)\) 평편을 통과할 때까지 수치적분을 수행한 후 평면을 수직으로 통과하는지 확인하면서 초기값 \(\mathbf{x}(t_0 )\) 를 보정해야 한다. 이 때 사용되는 방법이 미분보정(differential correction)이다.

 

 

편의상 식 (1)의 운동방정식을 다음과 같이 벡터 형식으로 표현한다.

 

\[ \dot{\mathbf{x}}(t)= \mathbf{f}(\mathbf{x}(t)) \tag{7} \]

 

그리고 초기값 \(\mathbf{x}(t_0 )\) 일 때의 궤적을. \(\mathbf{x}(t; \mathbf{x}(t_0 ))\) 로 표현한다. 그리고 약간의 섭동이 포함된 초기값 \(\mathbf{x}(t_0 )+ \delta \mathbf{x}(t_0 )\) 일 때 시간 \(t+\delta t\) 까지 전개한 궤적을 \(\mathbf{x}(t+ \delta t; \mathbf{x}(t_0 )+ \delta \mathbf{x}(t_0 ))\) 로 표현하자. 그러면 두 궤적 사이의 차이 \(\delta \mathbf{x}(t+ \delta t)\) 는 다음과 같다.

 

\[ \delta \mathbf{x}(t+ \delta t)=\mathbf{x}(t+ \delta t; \mathbf{x}(t_0 )+ \delta \mathbf{x}(t_0 ))- \mathbf{x}(t; \mathbf{x}(t_0 )) \tag{8} \]

 

시간 \(t_1+ \delta t_1\) 에서 \(\delta \mathbf{x}(t_1+\delta t_1)\) 을 계산하기 위해서 위 식을 테일러 시리즈를 1차항에서 절삭하면 다음과 같이 된다.

 

\[ \begin{align} \delta \mathbf{x}(t_1+ \delta t_1 ) &= \mathbf{x}(t_1; \mathbf{x}(t_0 ))+ \frac{ \partial \mathbf{x}(t_1, \mathbf{x}(t_0 ))}{\partial \mathbf{x}(t_0 ) } \delta \mathbf{x}(t_0 ) \tag{9} \\ \\ & \ \ \ + \frac{ \partial \mathbf{x}(t_1, \mathbf{x}(t_0 ))}{\partial t } \delta t_1 + HOT - \mathbf{x}(t_1; \mathbf{x}(t_0 )) \\ \\ & \approx \frac{ \partial \mathbf{x}(t_1, \mathbf{x}(t_0 ))}{\partial \mathbf{x}(t_0 ) } \delta \mathbf{x}(t_0 ) + \frac{ \partial \mathbf{x}(t_1, \mathbf{x}(t_0 ))}{\partial t } \delta t_1 \end{align} \]

 

여기서 \( \frac{ \partial \mathbf{x}(t_1, \mathbf{x}(t_0 ))}{\partial t } \)은 \(\dot{ \mathbf{x}}(t_1 )= \mathbf{f}(\mathbf{x}(t_1 ))\) 과 같고, \( \frac{ \partial \mathbf{x}(t_1, \mathbf{x}(t_0 ))}{\partial \mathbf{x}(t_0 ) } \) 은 \(t=t_1\) 에서 계산한 값으로서 \(\delta t_1=0\) 일 때의 상태천이행렬 \(\Phi (t_1, t_0 )\) 과 같다. 상태천이행렬 \(\Phi (t_1, t_0 )\) 은 식 (7)을 궤적 \(\mathbf{x}(t)\) 에서 선형화한 미분방정식의 상태천이행렬이다 (https://pasus.tistory.com/276).

 

\[ \delta \dot{\mathbf{x}}(t)= \left. \frac{ \partial \mathbf{f} }{ \partial \mathbf{x} } \right|_{\mathbf{x}(t)} \delta \mathbf{x}(t) \tag{10} \]

 

이제 초기값 식 (5)에서 출발한 궤적이 시간 \(t_1=\frac{T}{2} \)에서 식 (6)으로 주어지는 주기궤도의 조건식에 부합하는 값이 될 수 있도록 미분보정을 이용하여 초기값을 조정해 보자. 초기값 \(\mathbf{x}(t_0 )+ \delta \mathbf{x}(t_0 )\) 가 주어졌을 때 궤적이 \((x-z)\) 평편을 통과하는 시간 \(t_1+ \delta t_1\) 에서 상태변수는 다음과 같이 계산된다.

 

\[ \begin{align} & \mathbf{x}(t_1+ \delta t_1; \mathbf{x}(t_0 )+ \delta \mathbf{x}(t_0 )) \tag{11} \\ \\ & \ \ \ \ \ \ \ = \mathbf{x}(t_1; \mathbf{x}(t_0 ))+ \Phi (t_1, t_0 ) \delta \mathbf{x}(t_0 )+ \mathbf{f}( \mathbf{x}(t_1 )) \delta t_1 \end{align} \]

 

식 (11)이 일반적인 미분보정 또는 슈팅방법에서 사용하는 식과 다른 점이 있는데 바로 \(\delta t_1\) 이다. 일반적인 미분보정에서는 최종시간 \(t_f\) 가 고정된 값이다. 하지만 여기서는 초기값을 변경하면 \((x-z)\) 평편을 통과하는 시간이 달라지기 때문에 최종시간은 상수가 아니라 변수다.

 

 

식 (5)로 주어지는 초기값이 \((x-z)\) 평편을 수직으로 통과한다는 보장이 없기 때문에 \(\mathbf{x}(t_1; \mathbf{x}(t_0 ))\) 는 다음과 같이 주어질 것이다.

 

\[ \mathbf{x}(t_1; \mathbf{x}(t_0 ))= \begin{bmatrix} x_1 \\ 0 \\ z_1 \\ v_{x1} \\ v_{y1} \\ v_{z1} \end{bmatrix} \tag{12} \]

 

여기서 \(y_1=0\) 인 이유는 궤적이 \((x-z)\) 평편을 통과하는 시간 \(t_1\) 에서의 값이기 때문이다. \((x-z)\) 평편을 수직으로 통과하려면 식 (6)의 주기궤도 조건식과 같이 \(v_{x1}=v_{z1}=0\) 으로 만들어야 한다. 그래서 초기값을 \(\delta \mathbf{x}(t_0 )\) 만큼 변경하려고 한다. \(\delta \mathbf{x}(t_0 )\) 는 다음과 같다.

 

\[ \delta \mathbf{x}(t_0 )= \begin{bmatrix} \delta x_0 \\ 0 \\ \delta z_0 \\ 0 \\ \delta v_{y0} \\ 0 \end{bmatrix} \tag{13} \]

 

그러면 \((x-z)\) 평편을 통과하는 시간도 \(\delta t_1\) 만큼 바뀔 것이다. 여기서 목표는 식 (11)에서 \(\mathbf{x}(t_1 + \delta t_1; \mathbf{x}(t_0 )+ \delta \mathbf{x}(t_0 ))\) 값이 다음과 같이 되도록 \(\delta \mathbf{x}(t_0 )\) 를 계산하는 것이다.

 

\[ \mathbf{x}(t_1 + \delta t_1; \mathbf{x}(t_0 )+ \delta \mathbf{x}(t_0 ) ) = \begin{bmatrix} x^* \\ 0 \\ z^* \\ 0 \\ v^*_y \\ 0 \end{bmatrix} \tag{14} \]

 

여기서 \((x^*, \ z^*, \ v_y^*) \) 는 미지의 값으로서 어떤 특정한 값이 되도록 하는게 목표가 아니다. 목표는 \(v_x^*=v_z^*=0\) 으로 만드는 것이다.

식 (11)에 의하면 \(\delta \mathbf{x}(t_0 )\) 은 다음식으로 계산할 수 있다.

 

\[ \delta \mathbf{x}(t_1+ \delta t_1 )= \Phi (t_1, t_0 ) \delta \mathbf{x}(t_0 )+ \mathbf{f}( \mathbf{x}(t_1 )) \delta t_1 \tag{15} \]

여기서

\[ \delta \mathbf{x} (t_1+ \delta t_1 )= \begin{bmatrix} x^* \\ 0 \\ z^* \\ 0 \\ v^*_y \\ 0 \end{bmatrix}-\begin{bmatrix}x_1 \\ 0 \\ z_1 \\ v_{x1} \\ v_{y1} \\ v_{z1} \end{bmatrix} =\begin{bmatrix} \delta x_1 \\ 0 \\ \delta z_1 \\ \delta v_{x1} \\ \delta v_{y1} \\ \delta v_{z1} \end{bmatrix} \]

 

이다. 위 식을 풀어 쓰면 다음과 같다.

 

\[ \begin{align} 0 &= \Phi_{21} \delta x_0+ \Phi_{23} \delta z_0+ \Phi_{25} \delta v_{y0}+ v_{y1} \delta t_1 \tag{16} \\ \\ \delta v_{x1} &= \Phi_{41} \delta x_0+ \Phi_{43} \delta z_0+ \Phi_{45} \delta v_{y0}+ \dot{v}_{x1} \delta t_1 \\ \\ \delta v_{z1} &= \Phi_{61} \delta x_0+ \Phi_{63} \delta z_0+ \Phi_{65} \delta v_{y0}+ \dot{v}_{z1} \delta t_1 \end{align} \]

 

여기서 \(\Phi_{ij}\) 는 상태천이행렬 \(\Phi (t_1, t_0 )\) 의 성분이다. \(\delta x_1, \ \delta z_1, \ \delta v_{y1}\) 식을 제외한 이유는 보정할 필요가 없기 때문이다. 식 (16)에 의하면 식이 3개, 미지수는 4개이다. 따라서 미지수를 줄이기 위하여 \(x_0\) 는 고정시킨다. 즉 \(\delta x_0=0\) 로 놓는다, 그러면 식 (16)은 다음과 같이 된다.

 

\[ \begin{align} 0 &= \Phi_{23} \delta z_0+ \Phi_{25} \delta v_{y0}+ v_{y1} \delta t_1 \tag{17} \\ \\ \delta v_{x1} &= \Phi_{43} \delta z_0+ \Phi_{45} \delta v_{y0}+ \dot{v}_{x1} \delta t_1 \\ \\ \delta v_{z1} &= \Phi_{63} \delta z_0+ \Phi_{65} \delta v_{y0}+ \dot{v}_{z1} \delta t_1 \end{align} \]

 

식 (17)에서 \((\delta z_0, \ \delta v_{y0})\) 을 계산한 후 초기값을 다음식으로 업데이트한 다음,

 

\[ \begin{align} z_0 \ & \gets \ z_0+ \delta z_0 \tag{18} \\ \\ v_{y0} \ & \gets \ v_{y0}+ \delta v_{y0} \end{align} \]

 

위 프로세스를 반복하면 주기궤도를 만드는 초기값 \(\mathbf{x}(t_0 )\) 를 계산해 낼 수 있다.