강체 역학을 위한 FEM: 쿼터니언 기반 회전 처리와 관성 텐서

사면체 12자유도는 강체 구속 뒤 6이다. 간격 상한은 2/ω_max이고, 중간 자세는 SLERP가 구면 위에 남긴다.

강체 역학의 기본 방정식#

유체-구조 연성(FSI)이나 DEM에서 강체는 두 점 사이 거리만 지키면 된다. 변형은 없다. 그 거리를 수치로 지키려면 회전이 먼저 무너지지 않아야 한다.

입자 $i$ 의 속도는 질량중심 병진과 회전의 합이다.

$\mathbf{v}_i = \mathbf{v}_{cm} + \boldsymbol{\omega} \times \mathbf{r}_i$

$\mathbf{v}_{cm}$ 은 질량중심 속도, $\boldsymbol{\omega}$ 는 각속도, $\mathbf{r}_i$ 는 질량중심에서 입자까지의 변위다.

축 슬라이더와 ω를 바꾸면, 각속도 크기가 같아도 회전축이 바뀌고 한 바퀴 시간이 달라진다.

axis x 0.40axis y 1.00axis z 0.20ω 1.2

주황 화살표가 회전축. 사원수 q를 직접 곱해 매 프레임 큐브를 회전시킨다 — 짐벌락 없이 안정.

ω를 키우면 같은 축에서 자세가 더 빨리 돈다. 축을 몸통 대각으로 눕혀도 각이 한 곳에 잠기지 않는지 본다.

선형 모멘텀과 각 모멘텀#

전체 선형 모멘텀은 뉴턴 제2법칙으로부터

$\mathbf{P} = M \mathbf{v}_{cm}, \qquad \frac{d\mathbf{P}}{dt} = \mathbf{F}_{ext}$

여기서 $M = \sum_i m_i$ 는 전체 질량이다. 질량 중심의 가속도는

$M \dot{\mathbf{v}}_{cm} = \mathbf{F}_{ext}$

전체 각 모멘텀은

$\mathbf{L} = \mathbf{I} \boldsymbol{\omega}, \qquad \frac{d\mathbf{L}}{dt} = \boldsymbol{\tau}$

여기서 $\mathbf{I}$ 는 관성 텐서(inertia tensor), $\boldsymbol{\tau}$ 는 전체 토크(torque)다.

관성 텐서#

관성 텐서는 질량 중심을 기준으로

$I_{jk} = \sum_i m_i \left( |\mathbf{r}_i|^2 \delta_{jk} - r_{i,j} r_{i,k} \right)$

연속체로 확장하면

$\mathbf{I} = \int_V \rho \left( |\mathbf{r}|^2 \mathbf{E} - \mathbf{r} \otimes \mathbf{r} \right) dV$

여기서 $\mathbf{E}$ 는 단위 텐서이고, $\mathbf{r} \otimes \mathbf{r}$ 은 외적(outer product)이다.

물체가 회전하면 관성 텐서는 방위 행렬(orientation matrix) $\mathbf{R}$ 에 의해 변환된다.

$\mathbf{I}(t) = \mathbf{R}(t) \, \mathbf{I}_0 \, \mathbf{R}^T(t)$

여기서 $\mathbf{I}_0$ 는 물체 좌표계(body frame)에서의 초기 관성 텐서로, 시뮬레이션 시작 시 한 번만 계산하면 된다.

쿼터니언을 이용한 회전 표현#

오일러 각(Euler angles)으로 회전을 표현하면 짐벌 락(gimbal lock) 문제가 발생한다. 이를 피하기 위해 **쿼터니언(quaternion)**을 사용한다.

쿼터니언 $\mathbf{q} = (q_0, q_1, q_2, q_3)$ 는 스칼라 부분과 벡터 부분으로 구성된다.

$\mathbf{q} = q_0 + q_1 \mathbf{i} + q_2 \mathbf{j} + q_3 \mathbf{k}$

단위 축 $\hat{\mathbf{n}}$ 에 대해 각도 $\theta$ 만큼의 회전은

$\mathbf{q} = \cos\frac{\theta}{2} + \sin\frac{\theta}{2}\left(n_x \mathbf{i} + n_y \mathbf{j} + n_z \mathbf{k}\right)$

두 쿼터니언 $\mathbf{q}_1$ , $\mathbf{q}_2$ 의 합성 회전은 곱으로 표현된다.

$\mathbf{q}_{12} = \mathbf{q}_2 \mathbf{q}_1$

이를 이용하면 방위 행렬 $\mathbf{R}$ 은 쿼터니언 성분으로 명시적으로 표현된다.

\mathbf{R} = \begin{pmatrix} 1 - 2(q_2^2+q_3^2) & 2(q_1 q_2 - q_0 q_3) & 2(q_1 q_3 + q_0 q_2) \\ 2(q_1 q_2 + q_0 q_3) & 1 - 2(q_1^2+q_3^2) & 2(q_2 q_3 - q_0 q_1) \\ 2(q_1 q_3 - q_0 q_2) & 2(q_2 q_3 + q_0 q_1) & 1 - 2(q_1^2+q_2^2) \end{pmatrix}

쿼터니언의 시간 미분은 각속도 $\boldsymbol{\omega}$ 와 연결된다.

$\dot{\mathbf{q}} = \frac{1}{2} \begin{pmatrix} 0 \\ \boldsymbol{\omega} \end{pmatrix} \otimes \mathbf{q}$

$I_{jk}$ 의 대각은 $m(|\mathbf{r}|^2-r_j^2)$ , 비대각은 $-m r_j r_k$ 다. 질량이 $x$ 축 위에만 있으면 $r_y=r_z=0$ 이라 $I_{xx}=0$ , $I_{yy}=I_{zz}=m x^2$ 다. 그 축으로는 돌릴 팔 길이가 없다. $\mathbf{I}_0$ 는 이 합을 물체 좌표에 한 번 저장한 값이다. 월드에서는 $\mathbf{R}\mathbf{I}_0\mathbf{R}^T$ 만 다시 만들면 된다. 물체 좌표 텐서에 월드 각속도를 그대로 곱하면 토크의 축이 빗나간다. 좌표를 옮길 때의 합동은 구성 텐서의 좌표변환에서 요소 강성에도 같은 꼴로 나온다.

$\theta=0$ 이면 $\cos(\theta/2)=1$ , $\sin(\theta/2)=0$ 이라 $\mathbf{q}=(1,0,0,0)$ 이다. 코드가 넣는 초기 쿼터니언과 같다. $\theta=\pi$ 이면 $\cos(\pi/2)=0$ , $\sin(\pi/2)=1$ 이라 스칼라가 죽고 벡터부만 축 $\hat{\mathbf{n}}$ 이다. $\theta=2\pi$ 이면 $\cos\pi=-1$ 이라 $\mathbf{q}=(-1,0,0,0)$ 이다. 회전은 항등과 같고 부호만 반대다. 한 바퀴를 쿼터니언은 두 번 덮는다. $\mathbf{R}$ 의 $(1,1)$ 성분은 $1-2(q_2^2+q_3^2)$ 라, 이 항등값에서 1이 되고 비대각은 0이다. $\dot{\mathbf{q}}$ 앞의 $1/2$ 는 각도 공식의 $\theta/2$ 와 같은 인자다. 그대로 더하면 $|\mathbf{q}|$ 가 1에서 밀린다. 선형 사면체는 절점 4, 절점당 변위 3이라 12자유도다. 강체 구속은 병진 3과 회전 3, 합 6으로 줄인다. $\sum_i L_i=1$ 은 좌표가 남을 조건이지, 그 6개를 혼자 만들지는 않는다.

FEM 이산화와 사면체 형상 함수#

FEM에서 강체의 변위장을 이산화할 때, 사면체 요소의 형상 함수(shape function)는

$N_i(\xi, \eta, \zeta) = L_i, \quad i = 1, 2, 3, 4$

여기서 $L_i$ 는 무차원 체적 좌표(barycentric coordinates)이며 $\sum_i L_i = 1$ 조건을 만족한다.

선형 사면체 요소에서 변위 보간은

$\mathbf{u}(\mathbf{x}) = \sum_{i=1}^{4} N_i(\mathbf{x}) \, \mathbf{u}_i$

강체 구속 조건을 적용하면 자유도는 6개(선형 3 + 회전 3)로 줄어든다.

Python 구현 예시#

import numpy as np
 
class RigidBody:
    def __init__(self, mass, I_body):
        self.M = mass
        self.I0 = I_body          # body frame 관성 텐서 (3x3)
        self.x_cm = np.zeros(3)   # 질량 중심 위치
        self.v_cm = np.zeros(3)   # 질량 중심 속도
        self.q = np.array([1., 0., 0., 0.])  # 쿼터니언 (w, x, y, z)
        self.omega = np.zeros(3)  # 각속도
 
    def quaternion_to_rotation(self):
        w, x, y, z = self.q
        return np.array([
            [1-2*(y**2+z**2),  2*(x*y - w*z),  2*(x*z + w*y)],
            [2*(x*y + w*z),  1-2*(x**2+z**2),  2*(y*z - w*x)],
            [2*(x*z - w*y),    2*(y*z + w*x),  1-2*(x**2+y**2)]
        ])
 
    def inertia_world(self):
        R = self.quaternion_to_rotation()
        return R @ self.I0 @ R.T
 
    def step(self, F_ext, tau_ext, dt):
        # 선형 운동 적분
        a_cm = F_ext / self.M
        self.v_cm += a_cm * dt
        self.x_cm += self.v_cm * dt
 
        # 회전 운동 적분
        I_world = self.inertia_world()
        alpha = np.linalg.solve(I_world, tau_ext - np.cross(self.omega, I_world @ self.omega))
        self.omega += alpha * dt
 
        # 쿼터니언 업데이트
        omega_q = np.concatenate([[0.], self.omega])
        q_dot = 0.5 * quaternion_multiply(omega_q, self.q)
        self.q += q_dot * dt
        self.q /= np.linalg.norm(self.q)  # 단위 쿼터니언 정규화
 
def quaternion_multiply(q1, q2):
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    return np.array([
        w1*w2 - x1*x2 - y1*y2 - z1*z2,
        w1*x2 + x1*w2 + y1*z2 - z1*y2,
        w1*y2 - x1*z2 + y1*w2 + z1*x2,
        w1*z2 + x1*y2 - y1*x2 + z1*w2
    ])

정규화를 빼면 R이 회전이 아니게 된다#

매 스텝 $|\mathbf{q}|$ 를 1로 나누지 않으면 $\mathbf{R}$ 의 열이 단위 길이를 잃는다. 코드의 q /= norm(q)가 그 나눗셈이다. 토크를 월드에서 구하고 $\mathbf{I}_0$ 를 그대로 쓰면 $\boldsymbol{\alpha}$ 가 다른 축을 가리킨다. 월드 텐서는 $\mathbf{I}_{world} = \mathbf{R} \mathbf{I}_0 \mathbf{R}^T$ 로 바꾼 뒤에야 $\mathbf{L}=\mathbf{I}\boldsymbol{\omega}$ 와 맞는다. 명시 적분의 상한은 $\Delta t < 2 / \omega_{max}$ 다. $\omega_{max}$ 가 크면 간격을 먼저 줄인다.

두 자세 사이는 곧은 선으로 잇지 않는다. $t=0$ 에서 $\mathbf{q}_0$ , $t=1$ 에서 $\mathbf{q}_1$ 인 식이 SLERP다.

$\mathbf{q}(t) = \mathbf{q}_0 \left(\mathbf{q}_0^{-1} \mathbf{q}_1\right)^t$

선형 보간은 그 구면 안으로 꺼진다. 회전 자유도를 요소 변형과 섞어 세면, 얇은 쉘의 전단 잠김처럼 구속이 운동을 삼킨다. 강체라면 자유도는 12가 아니라 6이다.