Inertia Tensor

Motivation

최근 어떤 제어 도중 transition function을 만들기 위해 자유공간에서 병진운동뿐만 아니라 회전을 포함한 역학을 사용할 일이 있었다. 이때 회전을 다루는 것은 생각보다 까다로워 다양한 내용을 공부하게 되었다. 이를 정리한다.

Rotation

회전은 근본적으로 어렵다. 3차원에서의 회전이 벡터공간을 이루지 않기 때문이다. 일반적으로 자연스럽게 느껴지거나 직관적으로 알 수 있는 다양한 수학적 성질들은 대부분 벡터 공간의 성질인 경우가 많다. 그러나 3차원 회전은 벡터공간이 아니므로 이러한 성질들이 적용되지 않아서 다루기 어렵다.

조금 더 구체적으로, 3차원 회전은 수학적으로 $SO(3)$ 군으로 나타낼 수 있다. 이 군의 성질은 다양한데, 이중 가장 아쉬운 것은 교환법칙이 성립하지 않는다는 것이다. 회전에는 스칼라곱 등의 다양한 연산을 잘 정의할 수 있다. 그러므로 교환법칙만 성립하면 바로 벡터공간이 되어서 벡터공간의 편리한 성질들을 이용할 수 있다. 그러나 $SO(3)$ 군은 교환법칙이 성립하지 않아서 어떠한 방법으로도 회전을 벡터와 같은 방식으로 편리하게 다룰 수 없다.

물론 오일러 각이나 쿼터니언은 벡터처럼 표현되기도 한다. 그러나 이는 단순히 표현일 뿐이며 벡터와 같이 덧셈의 교환법칙 등 유용한 성질은 갖지 않는다.

이로부터 회전은 수학적 표기와 그 물리적 의미가 상당히 무관하게 된다. 예컨대 어떤 방식으로 회전 $A$를 표현했더니 $[1,2,3]$이었고, 같은 방식으로 회전 $B$를 표현했더니 $[4,5,6]$이었다고 하자. 이때 어떤 표헌법을 택하든 $A+B$는 일반적으로 $[5,7,9]$가 아니다. 물론 교환법칙이 성립하지 않기 때문에 애초에 일반적으로 $A+B\neq B+A$이다.

더욱이 점질량으로 모델링할 수 있는 병진운동과 다르게 회전운동의 경우 물체의 질량 분포에 크게 영향을 받는다. 즉, 애초에 물체의 질량 분포를 무시하고 회전을 수학적으로 모델링하는 것 자체가 불가능하다. 예컨대 3차원상에서 회전하는 물체는 어떠한 외력이나 공기저항, 마찰 없이도 스스로 그 회전축이 달라지는 등 다양한 현상을 보인다.

대표적으로 테니스 라켓 정리로도 알려진 중간축 정리(Intermediate Axis Theorem)는 물체가 중간 크기의 관성 모먼트를 가진 축을 따라 회전할 때 불안정해서 다른 축을 기준으로 같이 회전하게 되는 현상을 설명하는데, 이 현상은 어떠한 외력이 없더라도 발생한다.

Modeling

이러한 이유로 임의의 물체를 고정된 개수의 서로 제약된 점들로 표현하는 것을 시도했다. 왜냐하면 이렇게 함으로써 라그랑주 역학을 통하여 회전이 없는 경우와 동일한 방법으로 회전을 다룰 수 있기 때문이다. 물론 라그랑주 역학은 회전이 포함된 시스템도 다를 수 있지만 구속된 입자계로 다루면 라그랑지언을 일반화할 수 있으므로 더 편리하고 논리적 오류가 발생할 가능성도 낮다.

Inertia Tensor

이를 위해 먼저 회전을 결정하는 간단한 물리량이 있는지를 조사했다. 만약 이런 물리량이 있다면 그 물리량을 동일하게 갖는 서로 구속된 점질량들로 물체를 모델링하면 되기 때문이다. 다중 진자와 같은 링크 구조는 이러한 방식으로 해석되는 대표적인 예시다. 반대로 이런 물리량이 없다면 그것은 이러한 방식의 모델링이 불가능하다는 의미로, 유체역학의 대부분의 문제가 이런 경우에 해당한다.

물론 회전에 대해서는 관성 모먼트가 알려져있는 물리량이다. 그러나 관성 모먼트는 특정한 축에 대해 정의되는 실수로써 이것만으로는 물체의 회전이 결정되지 않는다는 것은 자명하다.

그래서 위키백과와 ChatGPT등을 조사해본 결과 물체의 회전 특성은 관성 텐서에 의해 완전히 결정된다는 것을 알게 되었다. 관성 텐서는 다음과 같이 정의되는 양이다.

$$I_{ij} = \sum_k m_k (r_k^2 \delta_{ij} - r_{ki} r_{kj})$$

이때 $m_k$는 물체를 이루는 $k$번째 입자의 질량이고, $r_k$는 $k$번째 입자의 위치 벡터다. $\delta_{ij}$는 크로네커 델타다.

관성 텐서는 추상적으로는 모든 축에 대한 관성 모먼트를 하나의 행렬로 표현한 것으로 볼 수 있으며, 다르게는 다음과 같이 각속도와 각운동량 사이의 관계를 결정하는 물리량으로도 생각할 수 있다.

$$L = I\omega$$

Normalization of Inertia Tensor

다만 관성 텐서 역시 물체의 위치(또는 좌표계 원점의 위치)와 물체의 회전(또는 좌표계의 회전)에 따라 그 값이 달라진다. 그러므로 이는 물체에 대한 내재적인 값이 아니다. 이러한 문제를 해결하기 위해 관성 텐서의 좌표계를 물체에 대하여 내재적인 값으로 정의할 필요가 있다. 먼저 물체의 형상에 대한 가장 기본적인 내재적 위치는 그 질량중심이다. 따라서 관성 텐서의 좌표계에서 물체의 질량중심을 원점으로 정의할 수 있다.

다음으로 물체의 방향을 결정해야 한다. 이때 다음과 같이 접근했다.

먼저 임의의 유한차원 실수 대칭 행렬은 스펙트럼 정리(Spectral Theorem)에 의해 항상 고유 벡터들로 이루어진 행렬과 고유값으로 이루어진 대각 행렬로 직교 대각화 가능(Orthogonal Diagonalizable)하며 그 고유값은 실수다. 즉, 그러한 대칭 행렬 $A$를 다음과 같이 고유 벡터로 이루어진 행렬 $Q$와 고유값으로 이루어진 대각 행렬 $\Lambda$로 분해할 수 있다.

$$A = Q\Lambda Q^T$$

행렬에 따라 고유값이 중복될 수 있다. 그러나 그런 경우에도 그람-슈미트 과정을 통해 직교 고유벡터 기저를 얻을 수 있다.

관성 텐서는 정의에 따라 대칭 행렬이므로 이러한 분해가 가능하다. 또한 고유 벡터들은 정의에 따라 $Iv = \lambda v$를 만족한다. 이때 $v$를 물체의 각속도라고 하면 $Iv$는 물체의 각운동량이 된다. 따라서 관성 텐서의 고유벡터는 그 축을 기준으로 물체를 회전시키면 각속도와 각운동량이 평행하게 되는 축을 의미한다. 이러한 축을 물체의 주축(principal axis)이라고 하며 그때의 고유벡터는 그 축을 기준으로 한 관성 모멘트가 된다.

또한 고유 벡터의 크기는 임의로 1로 설정할 수 있으므로 그렇게 하면 $Q$는 회전 행렬이 된다. 그러므로 임의의 관성 텐서에 대하여 그 고유 벡터로 이루어진 회전 행렬을 찾아서 물체에 적용하면 (또는 물체의 좌표계를 회전시키면) 관성 텐서를 대각 행렬로 만들 수 있다. 이렇게 하면 각 기저가 물체의 주축이 되며 각 기저의 고유값은 그 주축을 기준으로 한 관성 모멘트가 된다.

마지막으로 고유 벡터의 순서를 고유값의 크기 순으로 정렬해주기만 하면 관성 텐서는 완전히 물체의 내재적인 값이 된다. 좌표계의 대칭성(우수좌표계인지 좌수좌표계인지)는 이에 무관하므로 임의로 설정할 수 있다.

정리하자면 물체의 질량중심을 원점으로 이동시키고 물체의 주축이 그 크기 순으로 기저가 되도록 방향을 정렬하면 물체의 내재적인 관성 텐서를 얻을 수 있다.

Particle System

다음으로 이와 동일한 관성 텐서를 갖는 입자계를 만들어야 한다.

이를 위해 처음에는 자유도 분석에 기반한 입자계를 만들어보려 했으나 결론은 실패했다. 간단히 설명하자면 다음과 같다. 질량을 가진 점은 개당 4 자유도를 갖는다. 이로부터 축별 질량 중심 제약조건 3개와 관성 텐서의 상삼각 성분 10개, 질량 합에 대한 제약조건을 합하면 총 14개의 제약조건이 만들어진다. 이로부터 최소 4개의 점이 있어야 함을 알 수 있다. 이로부터 적당히 $n>4$에 대해 다음 연립방정식을 풀면 된다.

$$\begin{aligned} &\sum_{i=1}^n m_i = \text{total\_mass}\\ &\frac{\sum_{i=1}^n m_i x_i}{\sum_{i=1}^n m_i} = 0,\quad \frac{\sum_{i=1}^n m_i y_i}{\sum_{i=1}^n m_i} = 0,\quad \frac{\sum_{i=1}^n m_i z_i}{\sum_{i=1}^n m_i} = 0, \\ &\sum_{i=1}^n m_i (y_i^2 + z_i^2) = I_{xx},\quad \sum_{i=1}^n m_i (x_i^2 + z_i^2) = I_{yy},\quad \sum_{i=1}^n m_i (x_i^2 + y_i^2) = I_{zz}, \\ &\sum_{i=1}^n m_i x_i y_i = I_{xy},\quad \sum_{i=1}^n m_i y_i z_i = I_{yz},\quad \sum_{i=1}^n m_i x_i z_i = I_{xz}. \end{aligned}$$

문제는 이러한 연립방정식을 푸는 것이 쉽지 않다는 것이다. 물론 수치적인 방법을 동원하면 쉽게 해결할 수 있지만 더 간단한 방법이 있다.

물체의 질량중심이 원점에 있음에 착안, 질량이 $m$인 점 6 개가 $x$축 상에서 $+d_x, -d_x$에, $y$축 상에서 $+d_y, -d_y$에, $z$축 상에서 $+d_z, -d_z$에 위치하도록 입자계를 구성한다. 이러한 입자계는 원점이 질량중심이며 주축을 따라 정렬되었고 $m, d_x, d_y, d_z$ 4의 자유도를 갖는다. 이때 원점과 주축으로 정렬된 물체는 관성 텐서의 대각 성분과 질량을 합해서 총 4개의 자유도를 가지므로 이러한 입자계는 임의의 관성 텐서를 갖는 물체를 모델링할 수 있을 가능성이 높다. 물론 자유도 분석은 변수의 범위 등을 고려하지 않고 자유 변수의 개수만을 고려한 것이므로 이것이 임의의 질량과 관성 텐서를 모델링할 수 있는지는 별도로 증명해야 한다.

이러한 입자계는 다음과 같은 관성 텐서를 갖는다. 주축을 따라 정렬되었으므로 대각이 아닌 성분은 모두 0이다.

$$\begin{aligned} I_{xx} &= 2m(d_y^2 + d_z^2) \\ I_{yy} &= 2m(d_x^2 + d_z^2) \\ I_{zz} &= 2m(d_x^2 + d_y^2) \\ I_{xy} &= I_{yz} = I_{xz} = 0 \end{aligned}$$

이때 원래 물체의 질량을 $M$이라 하면 $m = M/6$이다. 이를 대입하면 다음을 얻는다.

$$\begin{aligned} I_{xx} &= \frac{M}{3}(d_y^2 + d_z^2) \\ I_{yy} &= \frac{M}{3}(d_x^2 + d_z^2) \\ I_{zz} &= \frac{M}{3}(d_x^2 + d_y^2) \end{aligned}$$

위 연립방정식을 풀기 위해 다음과 같이 이를 행렬 방정식으로 나타내어 풀었다.

$$\begin{bmatrix} I_{xx} \\ I_{yy} \\ I_{zz} \end{bmatrix} = \frac{M}{3} \begin{bmatrix} 0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \end{bmatrix} \begin{bmatrix} d_x^2 \\ d_y^2 \\ d_z^2 \end{bmatrix}\\ \therefore \begin{bmatrix} d_x^2 \\ d_y^2 \\ d_z^2 \end{bmatrix} = \frac{3}{M} \begin{bmatrix} 0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \end{bmatrix}^{-1} \begin{bmatrix} I_{xx} \\ I_{yy} \\ I_{zz} \end{bmatrix} = \frac{3}{2M} \begin{bmatrix} -1 & 1 & 1 \\ 1 & -1 & 1 \\ 1 & 1 & -1 \end{bmatrix} \begin{bmatrix} I_{xx} \\ I_{yy} \\ I_{zz} \end{bmatrix}$$

이때 좌변 벡터의 원소는 제곱수이므로 음수여서는 안 된다. 그런데 우변에서 만약 $I_{yy} + I_{zz} - I_{xx} < 0$이면 이는 음수가 될 수 있다. 따라서 그러한 상황이 발생하지 않음을 보여야 한다. 물론 이것은 수학적 모델이므로 허수 거리를 둘 수도 있지만, 그 경우 라그랑주 역학이 올바르게 동작할지 장담할 수 없다.

이는 다음과 같이 보일 수 있다. 먼저 단일 입자를 고려할 때 그 관성 텐서의 성분은 다음과 같이 주어진다.

$$I_{ij} = m (r^2 \delta_{ij} - r_i r_j)$$

알고자 하는 것은 대각 성분이므로 다음을 얻는다.

$$\begin{aligned} I_{xx} &= m (y^2 + z^2) \\ I_{yy} &= m (x^2 + z^2) \\ I_{zz} &= m (x^2 + y^2) \end{aligned}$$

이로부터 임의의 축 순서에 대하여

$$I_{jj}+I_{kk}-I_{ii} = m ((i^2 + k^2) + (i^2 + j^2) - (j^2 + k^2))=2mi^2 > 0$$

이때 여러 입자로 이루어진 물체의 관성 텐서는 각 입자의 관성 텐서의 합이고 부등식은 합에 대해 유지되므로 전체 관성 텐서에서도 이러한 성질이 성립한다. 따라서 위 방법으로 물체의 관성 텐서를 모델링할 수 있다.

아래는 이를 시각화한 것이다.

위 시뮬레이션은 화면이 작은 경우 제대로 보이지 않을 수 있으므로 이 링크를 클릭하여 전체화면으로 보는 것을 권장한다.

System

다음으로 이에 기반한 외력이 없는 상황에서의 물체의 회전을 라그랑주 역학으로 모델링하고자 시도했다. 처음에 이를 시도할 때에는 라그랑주 역학 부분은 간단하게 해결될 것이라 생각했는데 실제로 해 보니 그렇지 않았다.

먼저 물체의 일반화 좌표로 오일러 각을 사용했다. 만약 물체의 위치가 변할 수 있다면 물체의 회전 뿐만 아니라 좌표도 필요하다. 그러나 이 시스템은 질량중심이 곧 원점이며 외력이 없는 경우를 가정하기 때문에 물체의 위치는 변하지 않는다. 그러므로 물체의 회전만을 고려하면 된다.

다음으로 라그랑지언을 구해야 한다. 그런데 중력을 포함한 외력이 없으므로 퍼텐셜은 없다. 그러므로 라그랑지언은 곧 운동에너지가 된다.

이때 계의 운동에너지는 다음과 같이 주어진다.

$$T = \frac{1}{2} \sum_i m_i v_i^2 = \frac{1}{2} \sum_i m_i \dot{r}_i^2$$

이 식은 크게 복잡하지 않다. 그런데 여기서 $r_i$는 물체의 위치 벡터이므로 곧 물체의 기존 위치 $r_i^0$와 회전 행렬 $R$의 곱으로 주어진다. 회전 행렬은 다음과 같다.

$$R = \begin{bmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \cos\phi & 0 & \sin\phi \\ 0 & 1 & 0 \\ -\sin\phi & 0 & \cos\phi \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\psi & -\sin\psi \\ 0 & \sin\psi & \cos\psi \end{bmatrix}$$

이것을 모두 이용하여 라그랑지언을 구하면 다음과 같다. 식이 복잡하므로 sympy 라이브러리를 이용하여 계산했다. 질량을 포함한 비례 상수는 무시했다.

$$\begin{align*} L&\propto\\ &+l_{x}^{2} \left(\cos^{2}{\left(\theta{\left(t \right)} \right)} \left(\frac{d}{d t} \psi{\left(t \right)}\right)^{2} + \left(\frac{d}{d t} \theta{\left(t \right)}\right)^{2}\right) \\ &+ l_{y}^{2} \left(\sin^{2}{\left(\phi{\left(t \right)} \right)} \sin^{2}{\left(\theta{\left(t \right)} \right)} \left(\frac{d}{d t} \psi{\left(t \right)}\right)^{2} - \sin^{2}{\left(\phi{\left(t \right)} \right)} \left(\frac{d}{d t} \psi{\left(t \right)}\right)^{2} + \sin^{2}{\left(\phi{\left(t \right)} \right)} \left(\frac{d}{d t} \theta{\left(t \right)}\right)^{2} - 2 \sin{\left(\phi{\left(t \right)} \right)} \cos{\left(\phi{\left(t \right)} \right)} \cos{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \psi{\left(t \right)} \frac{d}{d t} \theta{\left(t \right)} - 2 \sin{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} \frac{d}{d t} \psi{\left(t \right)} + \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2} + \left(\frac{d}{d t} \psi{\left(t \right)}\right)^{2}\right) \\ &+ l_{z}^{2} \cdot \left(2 \sin{\left(\phi{\left(t \right)} \right)} \cos{\left(\phi{\left(t \right)} \right)} \cos{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \psi{\left(t \right)} \frac{d}{d t} \theta{\left(t \right)} - 2 \sin{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} \frac{d}{d t} \psi{\left(t \right)} - \cos^{2}{\left(\phi{\left(t \right)} \right)} \cos^{2}{\left(\theta{\left(t \right)} \right)} \left(\frac{d}{d t} \psi{\left(t \right)}\right)^{2} + \cos^{2}{\left(\phi{\left(t \right)} \right)} \left(\frac{d}{d t} \theta{\left(t \right)}\right)^{2} + \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2} + \left(\frac{d}{d t} \psi{\left(t \right)}\right)^{2}\right) \end{align*}$$

이제 일반화 좌표 $\theta, \phi, \psi$에 대해 다음 오일러-라그랑주 방정식을 연립하여 풀면 된다.

$$\frac{d}{dt} \left(\frac{\partial L}{\partial \dot{q}_i}\right) - \frac{\partial L}{\partial q_i} = 0$$

그런데 이것을 계산했더니 식을 정리하고 푸는 데 5분 이상이 소요되며 그 결과 역시 텍스트로 식을 저장한 것만 수 킬로바이트나 되었다. 그래서 블로그 포스트에 그것을 넣기는 어렵고 대신 다음 텍스트 파일을 첨부한다.

• 라그랑지안 계산 결과

이것은 다음과 같은 코드로 계산했다.

• 라그랑지안 및 시스템 방정식 계산 코드

이후 이것이 실제로 시스템을 모델링하는 데 사용될 수 있는지 알아보기 위해 간단한 시스템을 만들어 중간축 정리를 시험해보았다. 아래는 그 영상이다.

영상에서 후반으로 갈수록 약간씩 속도가 빨라지는데, 이것은 수치적인 오류로 인한 것으로 생각된다.

이 영상은 다음과 같은 코드로 만들어졌다.

• 시스템 시뮬레이션 코드

Conclusion

물체의 회전에 대해 공부했다. 적어도 모든 강체는 6개의 대칭적인 질량점의 조합으로 모델링 될 수 있음을 보였고, 이것을 시각화하여 직관과 일치하는지 확인했다. 또한 이러한 시스템이 주어졌을 때 라그랑주 역학을 사용하여 그 시스템 방정식을 얻어낼 수 있음을 보였다.