rigid body collision 1

물리 시뮬레이션 엔진에서의 기초적인 처리는 '충돌'이 아닐까 생각한다. 충돌을 처리하는 작업은 충돌하는 두 대상의 성질과 주변 환경의 작용에 따라서 쉬울수도, 굉장히 어려울수도 있다. 시뮬레이션에서는 모든 변수를 고려할 수 없기 때문에, 제한적인 상황을 가정하고 시뮬레이션 하는 경우가 대부분이다. 몇 번의 포스팅에 걸쳐 두 물체의 충돌에 대해서 설명해보려 한다.

 

문제를 간략히 하기위해 충돌하는 두 물체는 강체(Rigid body)라고 가정한다. 강체는 힘에 의해 형체가 변하지 않는 물체인데, 현실에선 존재하지 않지만 물리법칙을 설명할 때 다른 변수를 통제해서 식을 간단하게 한다. 물론 현실에서도 딱딱한 나무나 돌덩어리 등 변형정도가 무시할 수 있을만큼 적은 물체라면 강체라고 가정하여도 큰 무리가 없다.

 

물체의 충돌을 다루기 전에, 단일한 물체에 임의의 힘이 작용하는 경우에 대해서 먼저 이해해야 한다. 이번 포스팅은 이 주제를 다룰것이다.

 

물체에 힘이 작용한다는 것은 무엇을 의미하는가? 아이작 뉴턴은 운동에 관한 세 가지 법칙을 통해 이를 명쾌하게 설명했다.

 

1. 관성의 법칙: 물체의 질량중심은 외부 힘이 작용하지 않는 한 일정한 속도로 움직인다.

2. 가속도의 법칙: 물체의 운동량의 시간에 따른 변화율은 그 물체에 작용하는 알짜힘과 같다. 즉 [math]\mathbf{F} = \Delta \mathbf{P} = m\mathbf{a}[/math]

3. 작용 반작용의 법칙: 물체 A가 물체 B에 힘 [math]\mathbf{F}[/math]를 가하면, B는 A에게 크기가 같고 방향이 반대인 힘 [math]-\mathbf{F}[/math]를 동시에 가한다.

 

흔히 뉴턴의 운동 법칙이라 부르는 이 법칙에서 가장 주목할 부분은 제 2법칙이다. 그 유명한 [math]\mathbf{F} = m\mathbf{a}[/math]!!

 

이제 외부에서 힘이 가해졌을 때 물체에 일어나는 변화에 생각해보자. 외부에서 작용하는 힘은 물체를 '이동(translation)'시키거나 '회전(rotation)' 시킬 수 있다. 이 두 가지 변화를 계산한다면 힘의 작용을 시뮬레이션 할 수 있는것이다.

 

물체에 대한 힘의 작용 전에 입자의 이동에 관하여 살펴보자. 이동은 특정한 시간에서의 '위치(position)'로 표현할 수 있으며 시간 [math]t[/math]에서의 입자는 [math]x(t)[/math]로 표현한다. 입자를 이동하게 만드는 '속도(velocity)' 또한 시간에 대하여 표현되는데 [math]v(t)[/math] 로 나타내어진다. 시간 [math]t[/math]에서 속도는 위치의 시간에 대한 미분값이며 [math]\dot{x}(t)[/math]로 표현하기도 한다.

 

시간 변화량에 대한 입자의 이동은 오일러방법(Euler method)을 통해 해결할 수 있다. [math]\Delta t[/math] 후의 위치는 [math]x(t + \Delta t)[/math] 이며

[math!]x(t + \Delta t) = x(t) + \dot{x}(t) \Delta t[/math!]

와 같다. 즉, 시간 [math]t[/math]에서의 입자의 속도를 안다면 [math]\Delta t[/math]후의 위치를 계산할 수 있는것이다.

사족이지만, 오일러방법은 테일러급수에서 파생된것으로 [math]\Delta t[/math]의 크기가 커질수록 부정확한 값을 가진다. 당연히 0에 수렴하는 시간변화에 대한 위치를 계속하여 구할 수 있으면 좋겠지만, 현실적으로 불가능하므로 여러가지 방법을 통해 이를 보완하는 기법이 존재한다. 기회가 되면 테일러급수로부터 이를 보완하는 방법을 소개하겠다.

 

속도 또한 위치와 마찬가지로 오일러방법을 통해 구해지며 속도의 시간에 대한 미분값인 '가속도(acceleration)'가 사용된다. [math]a(t)[/math] 혹은 [math]\dot{v}(t)[/math]로 나타낸다. 이를 통해 [math]\Delta t[/math]후의 속도를 구하면 아래와 같다.

[math!]v(t + \Delta t) = v(t) + \dot{v}(t) \Delta t[/math!]

그렇다면 가속도는 어떻게 구한단 말인가? 처음 설명하길 외부에서 힘이 가해질 때 물체의 이동에 대하여 살펴보겠다고 했다. 즉, 뉴턴의 제 2법칙으로부터 가속도가 구해지는것이다.

[math!]a(t) = \frac{F}{m}[/math!]

이로부터 입자의 이동이 설명되었다. 하지만 입자와는 다르게 물체는 '회전'의 요소 또한 가지고 있으며 이동과 함께 이를 고려하여야 한다.

 

 

회전을 주관하는 속도는 '각속도(angular velocity)'라 부르며 [math]\omega (t)[/math]로 나타낸다. 각속도는 벡터로 그 방향은 회전의 축을 의미한다. 시간 [math]t[/math]에서 물체 위의 한 점 [math]s_i(t)[/math]에 작용하는 각속도는 회전반경의 접선방향으로 점을 회전시키는데, 이때 점의 속도는 물체의 '무게중심(center of mass)'으로부터 점까지의 거리 [math]r_i(t) = s_i(t) - CM[/math]와의 외적으로 표현된다.

[math!]\dot{ s_i}(t) = \omega(t) \times r_i(t)[/math!]

물체의 질량은 모든 점질량의 합이며, 무게중심은 물체 질량에 대한 각 점질량과 점의 위치의 곱의 합이다.

[math!]\begin{aligned} m &= \Sigma m_i\\ CM &= \frac{\Sigma m_i s_i}{m} \end{aligned}[/math!]

또한 각속도의 크기(magnitude)는 시간에 대한 회전각을 나타낸다. 즉, [math]|\omega(t)| = \theta (radian/sec)[/math]. 회전으로 인한 호의 길이가 [math]arc = r_i(t)\theta[/math] 이므로 [math]\omega(t) \times r_i(t)[/math]가 점의 속도를 나타낸다는 것을 이해할 수 있다.

 

이동과 같은 이치로 각속도를 통해 회전요소가 변화되고(이부분은 조금 후에 설명하겠다.) 각속도 또한 '각가속도(angular acceleration)'에 의해 변화된다. 마찬가지로 각가속도를 구하기 위해 회전에 작용하는 힘을 구해야 하는데 주어진 외력은 회전힘과 이동힘을 따로 구분하지 않았다. 물체 위 한 점 [math]s_i(t)[/math]에 작용하는 힘은 무게중심 방향으로 작용하는 힘과 이에 수직하는 접평면에 평행한 힘으로 분해될 수 있다.

[math!]f_i = f_i' + f_i^T[/math!]

[math]f_i^T[/math]가 회전을 일으키는 힘이며 이로부터 점에 작용하는 '돌림힘(torque)'이 계산된다.

[math!]\tau_i(t) = r_i(t) \times f_i^T[/math!]

그런데 물체에 하나의 힘만 작용하는것은 아니므로 모든 작용하는 힘으로부터 물체에 작용하는 돌림힘이 계산된다.

[math!]\tau(t) = \Sigma \tau_i(t) = \Sigma r_i(t) \times f_i^T[/math!]

돌림힘은 알았는데, 뉴턴의 제 2법칙이 어떻게 적용된단 말인가? [math]\mathbf{F} = m \mathbf{a}[/math] 보다 '운동량의 시간에 따른 변화율'에 주목해보면 답을 알 수 있다.

 

'선운동량(linear momentum)'은 모든 점질량과 속도의 곱으로, [math]\mathbf{P}(t) = \Sigma m_i v(t)[/math] 로 정의된다. 그런데 모든 점질량의 합은 물체의 질량을 나타내고 모든 점에서의 속도는 일정하므로

[math!]\mathbf{P}(t) = Mv(t)[/math!]

라고 나타낼 수 있으며, 이로부터

[math!]\mathbf{F}(t) = \dot{ \mathbf{P}(t) } = m \dot{ v(t) } = ma(t)[/math!]

임을 알 수 있는것이다.

 

'각운동량(angular momentum)'도 비슷한 맥락을 가진다. 차이점은 물체 위 한 점에서 각운동량은 무게중심으로부터의 거리[math]r_i(t)[/math]에 대한 운동량으로 표현되는것이다.

[math!]l_i(t) = r_i(t) \times m_i \left(\omega(t) \times r_i(t) \right) [/math!]

그리고 모든 점에 작용하는 각운동량이 물체의 각운동량을 정의한다.

[math!]L(t) = \Sigma l_i(t) = \Sigma r_i(t) \times m_i \left(\omega(t) \times r_i(t) \right)[/math!]

선운동량과 마찬가지로 돌림힘은 각운동량의 시간에 따른 변화율이다.

이를 쓰기에 앞서 각운동량을 조금 변형해보자.

[math!]\begin{equation}\begin{split} L(t) &= \Sigma r_i(t) \times m_i \left(\omega(t) \times r_i(t) \right)\\ &= \Sigma m_i r_i(t) \times \left(\omega(t) \times r_i(t) \right)\\ &= \Sigma -m_i r_i(t) \times r_i(t) \times \omega(t) \end{split}\end{equation}[/math!]

나름 깔끔하지만 만족할 수 없다. 여기에 두 벡터의 외적은 아래와 같이 메트릭스의 곱으로 표현될 수 있음을 적용해보자.

[math!]p \times q = \begin{bmatrix}0 & -p_z & p_y \\ p_z & 0 & -p_x \\ -p_y & p_x & 0 \end{bmatrix} \begin{bmatrix}q_x\\q_y\\q_z \end{bmatrix}[/math!]

즉 각운동량은 다음과 같이 다시 표현된다. ([math]r_{ix}(t)[/math] 등은 [math]r_{ix}[/math]로 간략히 썼습니다. 수식 자주쓰시는분들 존경스럽습니다..)

[math!]\begin{equation} \begin{split} L(t) &= \Sigma -m_i r_i(t) \times r_i(t) \times \omega(t) \\ &= \Sigma -m_i \begin{bmatrix} 0&-r_{iz} & r_{iy} \\ r_{iz} & 0 & -r_{ix} \\ -r_{iy} & r_{ix} & 0 \end{bmatrix} \begin{bmatrix} 0&-r_{iz} & r_{iy} \\ r_{iz} & 0 & -r_{ix} \\ -r_{iy} & r_{ix} & 0 \end{bmatrix} \omega(t) \\ &= \left( \Sigma m_i \begin{bmatrix}r_{iy}^2 + r_{iz}^2 & -r_{ix}r_{iy} & -r_{ix}r_{iz}\\ -r_{ix}r_{iy} & r_{iz}^2 + r_{ix}^2 & -r_{iy}r_{iz}\\ -r_{ix}r_{iz} & -r_{iy}r_{iz} & r_{ix}^2 + r_{iy}^2  \end{bmatrix} \right) \omega(t) \\ &= I(t) \omega(t) \end{split} \end{equation}[/math!]

선운동량이 [math]\mathbf{P}(t) = mv(t)[/math]임을 생각할 때, [math]I(t)[/math]는 [math]m[/math]와 비슷한 역할, 즉 '저항'의 역할을 하는 요소로 '관성텐서(inertia tensor)'라고 불린다. 선형힘과 같은 법칙이 돌림힘에 적용되며 관성텐서로 이를 표현할 수 있다.

[math!]\tau(t) = \dot{L}(t) = I(t)\dot{ \omega }(t) = I \alpha(t) [/math!]

 

다시 각가속도를 구하는 상황으로 돌아와 보자. 바로 위에서 장전된 사실로부터 각가속도를 구할 수 있다!!

[math!]\alpha(t) = I^{-1}(t) \tau(t)[/math!]

돌림힘은 [math]\tau(t) = \Sigma r_i(t) \times f_i[/math]로 부터 계산되므로 결국 초기에 작용하는 힘으로부터 각가속도를 계산할 수 있는것이다. 나아가 각가속도는 각속도를, 각속도는 회전요소를 변화시키게 되는것이다.

 

끝난것 같지만.. 아직 회전요소를 설명하지 않았다. 시간 [math]t[/math]에서의 회전요소는 3x3의 회전행렬 [math]\mathbf{R}(t)[/math]로 표현된다. 이제 오일러방법으로 회전요소만 갱신해주면 끝나는데..

[math!]\mathbf{R}(t + \Delta t) = \mathbf{R}(t) + \omega(t) \Delta t[/math!]

위 식은 성립할수가 없다. [math]\mathbf{R}(t)[/math]는 3x3행렬인 반면 [math]\omega(t)[/math]는 3차원 벡터인것이다. 이를 해결하기 위해 '회전사원수(rotation quaternion)'의 개념을 적용한다. (회전사원수는  rotate model using quaternion 에 설명했었습니다.) 회전요소는 오브젝트 공간의 물체를 월드공간의 회전상태로 만드는 요소이며 각속도 [math]\omega(t)[/math]에 의해 '더' 회전되어져 [math]\mathbf{R}(t + \Delta t)[/math]가 되는것임을 먼저 기억해두자.

회전사원수는 그와 동일한 회전을 정의하는 행렬로 변환될 수 있는것 처럼 역으로 회전행렬은 그와 동일한 회전을 정의하는 회전사원수로 변환될 수 있다. 즉, [math]\mathbf{R}(t)[/math]를 회전사원수 [math]q(t)[/math]로 변환하고, 각속도 [math]\omega(t)[/math]를 회전사원수 [math]q_{\omega}(t)[/math]로 변환하여 [math]\Delta t[/math]후의 회전행렬의 회전사원수를 계산할 수 있다.

[math!]q(t + \Delta t) = q_{\omega}(t)q(t)[/math!]

[math]q(t + \Delta t)[/math]를 다시 3x3 회전행렬로 변환하여 최종 회전행렬인 [math]\mathbf{R}(t + \Delta t)[/math]를 계산하게 된다.

 

이로써 외부 힘에 의한 물체의 변화 흐름을 모두 설명하였다. 이러한 무형의 힘은 중력이나 제다이의 포스(?!)에 의한 물체의 상태변화를 시뮬레이션 하기에 충분하다. 하지만 한가지 추가적으로 설명할게 남아있다.

 

선형 가속도의 계산을 위한 질량 [math]m[/math]과 달리, 각가속도의 계산에 필요한 관성텐서 [math]I(t)[/math]는 매 프레임마다 새롭게 계산되어야한다. 이는 굉장히 부담스러운 작업일 수 밖에 없는데, 이를 상수화 해보자.

[math!]\begin{equation}\begin{split} I(t) &= \Sigma m_i \begin{bmatrix}r_{iy}^2 + r_{iz}^2 & -r_{ix}r_{iy} & -r_{ix}r_{iz}\\ -r_{ix}r_{iy} & r_{iz}^2 + r_{ix}^2 & -r_{iy}r_{iz}\\ -r_{ix}r_{iz} & -r_{iy}r_{iz} & r_{ix}^2 + r_{iy}^2  \end{bmatrix} \\ &= \Sigma m_i \left( |r_i|^2 \mathbf{I} - \begin{bmatrix}r_{ix}^2 & r_{ix}r_{iy} & r_{ix}r_{iz}\\ r_{ix}r_{iy} & r_{iy}^2 & r_{iy}r_{iz}\\ r_{ix}r_{iz} & r_{iy}r_{iz} & r_{iz}^2 \end{bmatrix} \right) \\ &= \Sigma m_i \left( |r_i|^2 \mathbf{I} - r_i r_i^T \right) \\ &= \Sigma m_i \left( |\mathbf{R}(t) r_{i0}|^2 \mathbf{I} - (\mathbf{R}(t) r_{i0})(\mathbf{R}(t) r_{i0})^T \right) \\ &= \Sigma m_i \left( |r_{i0}|^2 \mathbf{I} - \mathbf{R}(t) r_{i0} r_{i0}^T \mathbf{R}^T(t) \right) \\ &= \Sigma m_i \left( \mathbf{R}(t)|r_{i0}|^2 \mathbf{I} \mathbf{R}^T(t) - \mathbf{R}(t) r_{i0} r_{i0}^T \mathbf{R}^T(t) \right) \\ &= \mathbf{R}(t) \left[ \Sigma m_i \left(|r_{i0}|^2 \mathbf{I} - r_{i0} r_{i0}^T \right) \right]  \mathbf{R}^T(t) \\ \end{split}\end{equation}[/math!]

두번째 식은 [math]|p|^2 = p_x^2 + p_y^2 + p_z^2[/math]로 부터 유도되고, 네번째 식은 [math]r_i(t) = \mathbf{R}(t)r_{i0}[/math]로 부터 유도된다([math]r_{i0}[/math]는 초기상태의 물체의 무게중심으로 부터 한 점 까지의 벡터이다). 다섯번째 식은 처음 물체는 강체라는 가정에 의해 무게중심으로부터 한 점 까지의 벡터의 크기는 변함없음과 전치행렬의 성질로부터 유도되었고, 여섯번째 식은 [math]\mathbf{R}(t)\mathbf{I}\mathbf{R}^T(t) = \mathbf{R}(t)\mathbf{R}^T(t) = \mathbf{I}[/math]로부터 유도되었다.

마지막 식을 보면 가운데 대괄호 안은 물체의 초기상태로부터 충분히 계산될 수 있음을 알 수 있다. 이로부터 관성텐서의 초기상수,

[math!]I_0 = \Sigma m_i \left(|r_{i0}|^2 \mathbf{I} - r_{i0} r_{i0}^T \right)[/math!]

를 미리 계산하여두어 시뮬레이션시의 계산량을 줄일 수 있다. 시뮬레이션의 매 프레임에서는 [math]I^{-1}(t)[/math]가 필요하므로, 즉,

[math!]I^{-1}(t) = ( \mathbf{R}(t) I_0 \mathbf{R}^T(t) )^{-1} = (\mathbf{R}^T(t))^{-1} I_{0}^{-1} \mathbf{R}^{-1}(t) = \mathbf{R}(t) I_{0}^{-1} \mathbf{R}^T(t)[/math!]

이므로 초기 계산시 [math]I_{0}^{-1}[/math]를 구하여 놓는다.

 

다음 포스팅을 통해 물체와 물체가 충돌할 때를 설명하도록 하겠다.

 

참고: 3d graphics for game programming