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

강체 역학의 기본 방정식#

유체-구조 연성(FSI) 해석이나 DEM(Discrete Element Method) 시뮬레이션에서 강체의 운동을 정확히 기술하는 것은 필수적이다. 강체는 변형이 없는 물체로, 임의의 두 점 사이 거리가 항상 일정하게 유지된다.

강체 내 임의의 입자 $i$의 속도는 다음과 같이 분해된다.

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

여기서 $\mathbf{v}_{cm}$은 질량 중심의 속도, $\boldsymbol{\omega}$는 각속도, $\mathbf{r}_i$는 질량 중심으로부터 입자 $i$까지의 변위 벡터이다.

선형 모멘텀과 각 모멘텀#

전체 선형 모멘텀은 뉴턴 제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}$은 쿼터니언 성분으로 명시적으로 표현된다.

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

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

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
    ])

실무 주의사항#

1. 쿼터니언 정규화 망각 수치 적분이 반복되면 $|\mathbf{q}| \neq 1$이 되어 회전 행렬이 비직교(non-orthogonal)가 된다. 매 타임스텝마다 q /= norm(q) 를 반드시 수행해야 한다.

2. 관성 텐서의 좌표계 혼동 $\mathbf{I}_0$는 물체 고정 좌표계(body frame)에서 정의된다. 월드 좌표계(world frame)에서 토크를 계산했다면 $\mathbf{I}_{world} = \mathbf{R} \mathbf{I}_0 \mathbf{R}^T$로 변환 후 적용해야 한다.

3. 타임스텝 선택 강체의 최고 고유진동수 $\omega_{max}$에 대해 $\Delta t < 2 / \omega_{max}$ 조건(CFL 유사 조건)을 만족해야 수치 안정성이 보장된다. 고속 회전 물체일수록 작은 $\Delta t$가 필요하다.

4. 쿼터니언 보간 (SLERP) 두 방위 사이를 보간할 때는 선형 보간(LERP) 대신 구면 선형 보간(SLERP)을 사용해야 등속 회전이 보장된다.

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

강체 역학과 FEM의 결합은 FSI, DEM, 다물체 동역학 시뮬레이션의 기반이 된다. 쿼터니언 기반 회전 표현을 올바르게 구현하면 짐벌 락 없이 임의의 3차원 회전을 안정적으로 시뮬레이션할 수 있다.