지구 관성 좌표계에서 회전하는 인공위성의 6-DOF equation of motion을 구해보자.
Requirements skills
- ECEF 좌표계
- ECI 좌표계
- Vector Operation
- Matrix Operation
- Kinematics
- State Space Representation
Kinematics
그림과 같이 ECEF 좌표계에서 움직이는 위성이 있다. 이때 ECI 좌표계에서 정의된 위성의 위치를 $\vec{r}_i=xi+yj+zk$라고 하자. ECI좌표계에서 정의된 값을 밑첨자i 라고 쓰겠다. 이때 Newton's 2nd Law는 관성 좌표계에서만 성립하므로, 위성의 운동방성식을 Inertia Frame인 ECI 좌표계를 기준으로 기술해야 한다. 지구의 자전에 의한 위성의 속도 변화를 포함한 운동 방정식은 아래와 같이 기술할 수 있다.
밑첨자 i는 ECI 좌표계, 밑첨자 e는 ECEF 좌표계를 의미한다.
$$\vec{v}_i = \frac{d}{dt}(\vec{r}_e)$$
$$=\dot{\vec{r}}_e+\vec{\omega}_E\times\vec{r}_e$$
$$\vec{a}_i = \frac{d^2}{dt^2}(\vec{r}_i)=\frac{d}{dt}(\vec{v}_i)$$
$$=\frac{d}{dt}(\dot{\vec{r}}_e+\vec{\omega}_E\times\vec{r}_e)$$
$$=(\ddot{\vec{r}}_e+\vec{\omega}_E\times\dot{\vec{r}}_e) + \frac{d}{dt}(\vec{\omega}_E)\times\vec{r}_e+\vec{\omega}_E\times\frac{d}{dt}(\vec{r}_e)$$
여기서 지구 자전 속도가 일정하다고 하면 $\frac{d}{dt}(\vec{\omega}_E)=0$ 이므로
$$=(\ddot{\vec{r}}_e+\vec{\omega}_E\times\dot{\vec{r}}_e) + \vec{\omega}_E\times(\dot{r}_e+\vec{\omega_E}\times\vec{r}_e)$$
$$=\ddot{\vec{r}}_e +2\vec{\omega}_E\times\dot{\vec{r}}_e+\vec{\omega}_E\times(\vec{\omega_E}\times\vec{r}_e)$$
그러므로 위성의 Kinematics는 아래와 같이 기술할 수 있다.
$$\vec{v}_i =\dot{\vec{r}}_e+\vec{\omega}_E\times\vec{r}_e $$
$$\vec{a}_i = \ddot{\vec{r}}_e +2\vec{\omega}_E\times\dot{\vec{r}}_e+\vec{\omega}_E\times(\vec{\omega_E}\times\vec{r}_e)=\Sigma{\frac{F_{ext}}{m}}$$
6-DOF Equation of Motion
위성의 운동 방정식에, 아래와 같이 정의된 벡터를 대입하자.
$$\vec{r}_e=xi+yj+zk$$
$$\dot{\vec{r}}_e=\dot{x}i+\dot{y}j+\dot{z}k$$
$$\ddot{\vec{r}}_e=\ddot{x}i+\ddot{y}j+\ddot{z}k$$
$$\vec{\omega}_E=\Omega_E k$$
$$\downarrow$$
$$2\vec{\omega}_E\times\dot{\vec{r}}_e=2(\Omega_E k)\times(\dot{x}i+\dot{y}j+\dot{z}k) = 2(\Omega_E\dot{x}j-\Omega_E\dot{y}i)$$
$$\vec{\omega}_E\times(\vec{\omega_E}\times\vec{r}_e)=\Omega_E k\times(\Omega_E k\times(xi+yj+zk))=\Omega_E k\times(\Omega_E xj-\Omega_E yi)$$
$$=-\Omega_E^2xi-\Omega_E^2yj$$
대입하여 정리하면 다음과 같이 6-DOF 방정식을 얻을 수 있다.
$$\therefore \ddot{x}i+\ddot{y}j+\ddot{z}k+2(\Omega_E\dot{x}j-\Omega_E\dot{y}i)-\Omega_E^2xi-\Omega_E^2yj=\Sigma{\frac{F_{ext}}{m}}$$
$$\begin{cases} \ddot{x} -2\Omega_E\dot{y}-\Omega_E^2x = F_x \\ \ddot{y} +2\Omega_E\dot{x}-\Omega_E^2y = F_y \\ \ddot{z}k =F_z \end{cases}$$
Gravity Force
위성에 작용하는 외력에는 중력과 추력이 있다. 중력은 아래와 같이 모델링할 수 있다.$($지구는 perfect shpere라고 가정한다.$)$
$$F_g = -\frac{\mu}{|r|^3}\vec{r} =-\frac{\mu}{(x^2+y^2+z^2)^{\frac{3}{2}}}(xi+yj+zk)$$
$(\mu = Gm_E = \textit{gravity coefficient})$
항력을 무시하고 추력을 사용하지 않는다고 가정하면 6DOF EOM은 아래와 같다.
$$\begin{cases} \ddot{x}-2\Omega_E\dot{y}-\Omega_E^2x=-\frac{\mu}{(x^2+y^2+z^2)^{\frac{3}{2}}}x \\ \ddot{y} +2\Omega_E\dot{x}-\Omega_E^2y=-\frac{\mu}{(x^2+y^2+z^2)^{\frac{3}{2}}}y \\ \ddot{z}k=-\frac{\mu}{(x^2+y^2+z^2)^{\frac{3}{2}}}z \end{cases}$$
State Space Representation
다음과 같이 위성의 State를 정의하자.
$$X=[x,y,z,\dot{x},\dot{y},\dot{z}]$$
이는 ECEF 좌표계를 기준으로 정의된 변수 이다.
$$\dot{X}=[\dot{x},\dot{y},\dot{z},\ddot{x},\ddot{y},\ddot{z}]$$
이때 위 연립 미분 방정식은 아래와 같이 나타낼 수 있다.
$$\begin{bmatrix}\dot{x} \\ \dot{y} \\ \dot{z} \\ \ddot{x} \\ \ddot{y} \\ \ddot{z} \end{bmatrix}=\begin{bmatrix} 0&0&0&1&0&0 \\ 0&0&0&0&1&0 \\ 0&0&0&0&0&1 \\ \Omega_E^2+A&0&0&0&2\Omega_E&0 \\ 0&\Omega_E^2+A&0&-2\Omega_E&0&0 \\ 0&0&A&0&0&0 \end{bmatrix} \begin{bmatrix} x\\y\\z\\ \dot{x} \\ \dot{y} \\ \dot{z} \end{bmatrix}$$
$$A = -\frac{\mu}{(x^2+y^2+z^2)^{\frac{3}{2}}}$$
이 식은 x,y,z가 서로 couple되어 있는 비 선형 연립 미분방정식 이다. 그러므로 이 식을 풀기 위해 Linearlization을 하거나, 계산 가능하도록 알고리즘을 짜야 한다.
만약 위성의 궤도가 Circular Orbit이라면 $x^2+y^2+z^2 = r(0) = \textit{constant}$라 가정하고 초기값을 대입하여 해를 구할 수 있다.
또는 수치 적분을 통해 step-by-step으로 시뮬레이션 할 수 있다.
cf$)$ 사용한 가정은 아래와 같다.
ⓛ Earth is a Perfect Sphere
② No Drag
③ No Thrust
④ No Perturbation
Drag - Decay Time
200km 이상의 고도에서는 밀도가 매우 낮지만, 결코 무시할 수는 없다. 위성에도 공기저항에 의해 항력이 작용하기 때문에 시간이 지날수록 속도가 감소하고 고도가 낮아진다. 항력은 다음과 같이 모델링한다.
$$\vec{D} = \frac{1}{2}\rho |\vec{v}|\vec{v}SC_d$$
S : Sectional Area
$C_d$ : Drag coefficient
$\rho$ : density
단면적과 항력계수는 위성의 형상에 따라 정해지므로 상수로 생각한다. 그러나 밀도는 고도에 따라 달라지므로 상수가 아니다. 위성의 현재 고도에 대한 함수로 밀도를 정의해 주어야 한다.
ISA table을 참고하면 특정 고도에서 밀도값을 알 수 있다. 피팅을 통해 밀도를 모델링할 수 있다.
피팅된 밀도를 $\rho(r)$라고 하면 항력에 의해 생기는 가속도는 다음과 같다.
$$\vec{a}_d = -\frac{\vec{D}}{m} = \frac{1}{2}\rho(r)|\vec{v}|\vec{v}\frac{SC_d}{m}$$
$$=-\frac{1}{2}\rho(r)\sqrt{x^2+y^2+z^2}(\dot{x}i+\dot{y}j+\dot{z}k)C_b$$
$C_b = \frac{SC_d}{m}$ : ballistic coefficient
그러므로 항력을 EOM에 추가하면 아래와 같다.
$$\begin{bmatrix}\dot{x} \\ \dot{y} \\ \dot{z} \\ \ddot{x} \\ \ddot{y} \\ \ddot{z} \end{bmatrix}=\begin{bmatrix} 0&0&0&1&0&0 \\ 0&0&0&0&1&0 \\ 0&0&0&0&0&1 \\ \Omega_E^2+A&0&0&B&2\Omega_E&0 \\ 0&\Omega_E^2+A&0&-2\Omega_E&B&0 \\ 0&0&A&0&0&B \end{bmatrix} \begin{bmatrix} x\\y\\z\\ \dot{x} \\ \dot{y} \\ \dot{z} \end{bmatrix}$$
$$A = -\frac{\mu}{(x^2+y^2+z^2)^{\frac{3}{2}}}$$
$$B = -\frac{1}{2}\rho(r)\sqrt{x^2+y^2+z^2}C_b$$
