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

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

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

 

$$\begin{array}{r}\text{(1)}& & \ddot{x}-2\dot{y}-x=-\frac{(1-\mu )(x+\mu )}{r_1^3}-\frac{\mu (x+\mu -1)}{r_2^3}\\ & \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{array}$$

여기서

$$\begin{array}{r}\text{(2)}& & r_1=\sqrt{(x+\mu {)}^2+y^2+z^2}\\ & r_2=\sqrt{(x+\mu -1{)}^2+y^2+z^2}\end{array}$$

 

이다.

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

 

$$\begin{array}{r}\text{(3)}& & \dot{x}=\frac{dx}{dt}=\frac{dx}{-d\tau }=-x^{\prime },\dot{y}=\frac{dy}{dt}=\frac{-dY}{-d\tau }=Y^{\prime }\\ & \ddot{x}=\frac{d^{2}x}{d{t}^2}=\frac{d^{2}x}{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}{d{t}^2}=\frac{d^{2}z}{d{\tau }^2}=z^{″}\end{array}$$

 

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

 

$$\begin{array}{r}\text{(4)}& & x^{″}-2Y^{\prime }-x=-\frac{(1-\mu )(x+\mu )}{r_1^3}-\frac{\mu (x+\mu -1)}{r_2^3}\\ -& Y^{″}-2x^{\prime }+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{array}$$

 

식 (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)$ 평편에서 수직으로 출발해야 하므로 주기궤도를 위한 초기조건은 다음과 같아야 한다.

 

$$\begin{array}{}\text{(5)}& \mathbf{x}(t_0)=\left[\begin{array}{c}x(t_0)\\ y(t_0)\\ z(t_0)\\ \dot{x}(t_0)\\ \dot{y}(t_0)\\ \dot{z}(t_0)\end{array}\right]=\left[\begin{array}{c}{x}_0\\ 0\\ z_0\\ 0\\ v_{y0}\\ 0\end{array}\right]\end{array}$$

 

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

 

$$\begin{array}{}\text{(6)}& \mathbf{x}(T/2)=\left[\begin{array}{c}x(T/2)\\ 0\\ z(T/2)\\ 0\\ v_y(T/2)\\ 0\end{array}\right]\end{array}$$

 

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

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

 

 

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

 

$$\begin{array}{}\text{(7)}& \dot{\mathbf{x}}(t)=\mathbf{f}(\mathbf{x}(t))\end{array}$$

 

그리고 초기값 $\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)$ 는 다음과 같다.

 

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

 

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

 

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

 

여기서 $\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$ 일 때의 상태천이행렬 $\mathrm{\Phi }(t_1,t_0)$ 과 같다. 상태천이행렬 $\mathrm{\Phi }(t_1,t_0)$ 은 식 (7)을 궤적 $\mathbf{x}(t)$ 에서 선형화한 미분방정식의 상태천이행렬이다 (https://pasus.tistory.com/276).

 

$$\begin{array}{}\text{(10)}& \delta \dot{\mathbf{x}}(t)={\frac{\partial \mathbf{f}}{\partial \mathbf{x}}|}_{\mathbf{x}(t)}\delta \mathbf{x}(t)\end{array}$$

 

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

 

$$\begin{array}{r}\text{(11)}& & \mathbf{x}(t_1+\delta t_1;\mathbf{x}(t_0)+\delta \mathbf{x}(t_0))\\ & =\mathbf{x}(t_1;\mathbf{x}(t_0))+\mathrm{\Phi }(t_1,t_0)\delta \mathbf{x}(t_0)+\mathbf{f}(\mathbf{x}(t_1))\delta t_1\end{array}$$

 

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

 

 

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

 

$$\begin{array}{}\text{(12)}& \mathbf{x}(t_1;\mathbf{x}(t_0))=\left[\begin{array}{c}{x}_1\\ 0\\ z_1\\ v_{x1}\\ v_{y1}\\ v_{z1}\end{array}\right]\end{array}$$

 

여기서 $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)$ 는 다음과 같다.

 

$$\begin{array}{}\text{(13)}& \delta \mathbf{x}(t_0)=\left[\begin{array}{c}\delta x_0\\ 0\\ \delta z_0\\ 0\\ \delta v_{y0}\\ 0\end{array}\right]\end{array}$$

 

그러면 $(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)$ 를 계산하는 것이다.

 

$$\begin{array}{}\text{(14)}& \mathbf{x}(t_1+\delta t_1;\mathbf{x}(t_0)+\delta \mathbf{x}(t_0))=\left[\begin{array}{c}{x}^{\ast }\\ 0\\ z^{\ast }\\ 0\\ v_y^{\ast }\\ 0\end{array}\right]\end{array}$$

 

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

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

 

$$\begin{array}{}\text{(15)}& \delta \mathbf{x}(t_1+\delta t_1)=\mathrm{\Phi }(t_1,t_0)\delta \mathbf{x}(t_0)+\mathbf{f}(\mathbf{x}(t_1))\delta t_1\end{array}$$

여기서

$$\delta \mathbf{x}(t_1+\delta t_1)=\left[\begin{array}{c}{x}^{\ast }\\ 0\\ z^{\ast }\\ 0\\ v_y^{\ast }\\ 0\end{array}\right]-\left[\begin{array}{c}{x}_1\\ 0\\ z_1\\ v_{x1}\\ v_{y1}\\ v_{z1}\end{array}\right]=\left[\begin{array}{c}\delta x_1\\ 0\\ \delta z_1\\ \delta v_{x1}\\ \delta v_{y1}\\ \delta v_{z1}\end{array}\right]$$

 

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

 

$$\begin{array}{r}\text{(16)}& 0& ={\mathrm{\Phi }}_{21}\delta x_0+{\mathrm{\Phi }}_{23}\delta z_0+{\mathrm{\Phi }}_{25}\delta v_{y0}+v_{y1}\delta t_1\\ \delta v_{x1}& ={\mathrm{\Phi }}_{41}\delta x_0+{\mathrm{\Phi }}_{43}\delta z_0+{\mathrm{\Phi }}_{45}\delta v_{y0}+{\dot{v}}_{x1}\delta t_1\\ \\ \delta v_{z1}& ={\mathrm{\Phi }}_{61}\delta x_0+{\mathrm{\Phi }}_{63}\delta z_0+{\mathrm{\Phi }}_{65}\delta v_{y0}+{\dot{v}}_{z1}\delta t_1\end{array}$$

 

여기서 ${\mathrm{\Phi }}_{ij}$ 는 상태천이행렬 $\mathrm{\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{array}{r}\text{(17)}& 0& ={\mathrm{\Phi }}_{23}\delta z_0+{\mathrm{\Phi }}_{25}\delta v_{y0}+v_{y1}\delta t_1\\ \delta v_{x1}& ={\mathrm{\Phi }}_{43}\delta z_0+{\mathrm{\Phi }}_{45}\delta v_{y0}+{\dot{v}}_{x1}\delta t_1\\ \\ \delta v_{z1}& ={\mathrm{\Phi }}_{63}\delta z_0+{\mathrm{\Phi }}_{65}\delta v_{y0}+{\dot{v}}_{z1}\delta t_1\end{array}$$

 

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

 

$$\begin{array}{r}\text{(18)}& z_0& \leftarrow z_0+\delta z_0\\ v_{y0}& \leftarrow v_{y0}+\delta v_{y0}\end{array}$$

 

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