Simple MATLAB simulation using Euler rigid body equation and Quaternion kinematics

이 포스팅에는 강체의 회전을 표현하는데 널리 사용되고 있는 오일러 강체 방정식( rigid body dynamics equation)과 사원수 기구학 방정식 (quaternion kinematics equation) 이용하여 Matlab으로 간단한 시뮬레이션 코드를 작성하는 것에대해 작성한다.

아래 그림은 시뮬레이션 코드 작성을 위해 간단히 표현한 블록도이다. 강체 방정식은 인공위성이나 어떠한 회전하는 물체가 어느정도의 토크를 받냐에 따라 회전력을 계산할 수 있다. 그리고 quaternion kinematic를 이용하여 동체의 회전 각도를 표현한다.

Figure 1. Dynamics and Kinematics System Block Diagram

Figure 1의 중앙 상단 블럭은 Euler rigid body equation을 풀어주는 부분이며 식은 다음과 같다.

$$\dot{\omega}=J_B^{-1}\left[\tau_{ext}-\omega\times J_B\omega\right]$$

$$\tt \omega:angular\:rate\:wrt\:body\:coordinate\:system\:\in R^{3x1}[m/sec]\\\tau_{ext}:external\:torque\:wrt\:body\:coordinate system\:\in R^{3x1}[Nm]\\J_B:moment\:of\:inertial\:matrix\:\in R^{3x3}[kgm^2]$$

블럭의 출력은 각속도의 변화량이므로 이를 적절한 적분 방법(ex ode, rk ...)으로 적분하여  중앙 하단 블럭의 Quaternion kinematics에 전달한다. 

이 블럭의 출력은 적절한 적분 방법을 거쳐 Quaternion의 요구조건을 충족하기위해 정규화를 진행하여 최종적으로 시뮬레이션 루프를 이룬다.

임의의 토크를 동체의 x축에 가했을 때의 시뮬레이션 결과는 다음과 같다. 결과는 당연히 가해진 x축을 기점으로 회전함을 알 수 있으며 특히 토크를 계속 가하기 때문에 각속도가 계속 증가함에 따라 자세각을 통해 계속 회전하고 회전속도가 빨라짐을 알수있다.

Matlab Code

Reference -

[1] https://en.wikipedia.org/wiki/Euler%27s_equations_(rigid_body_dynamics)

[2] https://phys.libretexts.org/Bookshelves/Classical_Mechanics/Variational_Principles_in_Classical_Mechanics_(Cline)/13%3A_Rigid-body_Rotation/13.17%3A_Eulers_equations_of_motion_for_rigid-body_rotation