Authors: Miles Macklin, Matthias Müller, Nuttapong ChentanezPublished in: MiG '16: Proceedings of the 9th International Conference on Motion in Games (2016)- 해당 포스팅은 AI를 통해 번역한 내용입니다. 일부 내용이나 해석에 오역 또는 의역이 있을 수도 있습니다.- 개인적으로 공부하며 핵심 내용을 정리하고 해설한 포스팅입니다.- 본 저작물의 원문 저작권은 저자 및 출판사에 있습니다. 문제가 될 시 삭제될 수 있으며, 원문은 영리적 목적으로 복사 혹은 배포가 금지되어 있음을 알립니다.
1. Introduction
위치 기반 역학(Position-Based Dynamics) [Müller 등 2007]은 게임 및 대화형 응용 프로그램(interactive applications)에서 변형 가능한 물체(deformable bodies)의 실시간 시뮬레이션(real-time simulation)을 위해 널리 사용되는 방법입니다. 이 방법은 특히 그 단순성(simplicity)과 강건성(robustness) 덕분에 매력적이며, 최근에는 게임 분야를 넘어 영화 및 의료 시뮬레이션 분야에서도 인기를 얻고 있습니다.
위치 기반 역학(PBD)의 인기가 높아짐에 따라 그 한계점들이 더욱 문제가 되고 있습니다. 잘 알려진 한 가지 한계는 PBD의 동작이 시뮬레이션의 시간 단계(time step)와 반복 횟수(iteration count)에 의존적이라는 점입니다 [Bender 등 2014b]. 구체적으로, 반복 횟수(iteration count)가 증가하거나 시간 단계(time step)가 감소함에 따라 제약 조건(constraints)이 임의로 단단해지게(stiff) 됩니다. 이러한 매개변수의 결합(coupling of parameters)은 다양한 재질 유형이 포함된 장면을 생성할 때 특히 문제가 되는데, 예를 들어 부드러운 물체(soft bodies)가 거의 강체(rigid bodies)에 가까운 물체와 상호작용하는 경우입니다.
기존 PBD(Position-Based Dynamics)는 계산 효율성은 높지만, 물리적 강성($k$)이 시뮬레이션 환경에 따라 변하는 치명적인 한계가 있습니다.
가장 큰 문제는 제약 조건을 만족시키기 위한 반복 횟수(Iteration count)가 많아질수록 물체가 의도치 않게 더 단단해지는 '역설'이 발생한다는 점인데, 이는 고사양 PC와 저사양 모바일 기기 간의 시뮬레이션 결과가 불일치하는 원인이 됩니다. 또한 강성 계수 $k$가 물리적 단위가 아닌 단순 수치라 시간 단계(Time step, $\Delta t$) 변화에 민감하게 반응하여 프레임레이트가 떨어지면 물체의 물리성이 무너지는 현상이 나타납니다.
이러한 시나리오에서 한 물체의 강성을 얻기 위해 반복 횟수를 높이면 시뮬레이션 내의 다른 모든 물체의 동작이 의도치 않게 변할 수 있습니다. 이는 종종 강성 계수(stiffness coefficients)를 전역적으로 다시 조정해야 하게 만들며, 재사용 가능한 시뮬레이션 자산(simulation assets)을 만드는 것을 매우 어렵게 합니다. 반복 횟수 의존성은 단일 자산의 경우에도 문제가 되는데, 예를 들어 옷감 모델에서 인장(stretch) 및 굽힘 제약(bending constraints)의 상대적 강성을 설정할 때가 그렇습니다. 설상가상으로 반복 횟수의 효과는 비선형적(non-linear)이어서, 매개변수를 직관적으로 조정하거나 강성 값을 단순히 반복 횟수의 함수로 재조정(rescale)하기 어렵게 만듭니다.
최근 가상 현실(virtual-reality)의 재부상은 더 높은 정밀도(fidelity)와 물리적으로 더 대변적인 실시간 시뮬레이션의 필요성을 불러일으켰습니다. 동시에, 촉각 피드백 장치(haptic feedback devices)의 광범위한 사용은 정확한 힘 추정치(force estimates)를 제공할 수 있는 방법을 요구합니다. PBD는 잘 정의된 제약 힘(constraint force)의 개념이 없으며, 이 때문에 정확성보다 속도가 더 중요하고 시뮬레이션이 부차적인 효과인 응용 분야에 주로 국한되어 왔습니다.
본 논문에서 우리는 확장된 위치 기반 역학(extended position-based dynamics, XPBD) 알고리즘을 제시합니다. 우리의 방법은 잘 정의된 탄성 위치 에너지(elastic potential energy) 개념에 상응하는 새로운 제약 조건 공식을 도입함으로써 반복 횟수 및 시간 단계에 의존적인 강성 문제를 해결합니다. 우리는 PBD에 전체 라그랑주 승수(total Lagrange multiplier) 개념을 도입하는 암시적 시간 이산화(implicit time discretization)로부터 우리의 방법을 도출합니다. 이는 힘에 의존하는 효과와 장치를 구동하는 데 사용할 수 있는 제약 힘 추정치를 제공합니다.
요약하자면, 우리의 주요 기여는 다음과 같습니다:
- 잘 정의된 탄성 및 소산 에너지 전위(dissipation energy potentials)와 직접적으로 대응하도록 PBD 제약 조건을 확장했습니다.
- PBD에 전체 라그랑주 승수(total Lagrange multiplier) 개념을 도입하여 시간 단계와 반복 횟수에 독립적인 방식으로 제약 조건을 해결할 수 있게 했습니다.
- 비선형 뉴턴 솔버(non-linear Newton solver) 기반의 기준 암시적 시간 단계화 기법(implicit time stepping scheme)과 비교하여 우리 알고리즘을 검증했습니다.
3 Background
우리는 위치 기반 역학(position-based dynamics) 알고리즘의 핵심을 간략하게 설명하며, 더 자세한 내용은 Bender 등 [2014b]의 서베이 논문을 참조하시기 바랍니다. PBD는 스토머-벌렛(Størmer-Verlet) 방법을 사용한 반암시적 적분(semi-implicit integration) 단계와 그 뒤를 잇는 다수의 제약 조건 투영(constraint projection) 단계로 생각할 수 있습니다.
- 반암시적 적분(Semi-implicit Integration)이란?
보통 물리 엔진은 '힘(Force) → 가속도 → 속도 → 위치' 순으로 계산합니다. 하지만 PBD는 스토머-벌렛(Størmer-Verlet) 방식을 빌려와, 가속도를 직접 다루기보다 이전 위치와 현재 위치의 차이를 통해 다음 위치를 예측합니다. 이 방식은 계산이 빠르면서도 에너지가 폭발하지 않고 안정적으로 유지되는 장점이 있습니다.
- 제약 조건 투영(Constraint Projection)의 직관적 이해
PBD의 독특한 점은 위치를 먼저 예측해 본 뒤, "어? 옷감이 너무 늘어났네?" 혹은 "공이 바닥을 뚫고 들어갔네?" 같은 제약 조건(Constraints) 위반을 발견하면, 그 자리에서 위치를 직접 옮겨버린다는 것입니다. 예시 : 실을 잡아당길 때 힘의 평형을 계산하는 것이 아니라, 실의 원래 길이보다 길어지면 강제로 점들을 원래 길이 위치로 끌어다 놓는 방식입니다.
- 왜 '다수의' 투영 단계가 필요한가?
한 번의 수정으로 모든 제약 조건을 만족시키기는 어렵습니다. 예를 들어, 한쪽 점을 옮기면 연결된 다른 쪽 실이 또 늘어날 수 있죠. 그래서 여러 번 반복(Iteration)하며 조금씩 모든 점을 제약 조건에 맞게 수렴(Convergence)시켜 나가는 과정이 필요합니다.
투영 단계는 질량 가중치 보정(mass weighted corrections)을 동반한 각 제약 함수(constraint function)의 국소적 선형화(local linearizations)를 사용하여 수행됩니다. PBD 제약 조건 솔버의 주요 단계는 다음과 같이 주어지는 제약 조건별 위치 변화량(position delta)의 계산입니다.
- (1) \Delta \mathbf{x} = kj sj \mathbf{M}-1 \nabla Cj (\mathbf{x}i) \quad
여기서 하첨자 $i$는 반복 인덱스(iteration index)를, $j$는 제약 조건 인덱스(constraint index)를 나타내며, $k \in [0, 1]$는 각 제약 조건 보정치에 단순한 곱수로 적용되는 제약 강성(constraint stiffness)입니다. 스케일링 인자 $s$는 제약 함수의 단일 뉴턴 단계(Newton step)로부터 유도된 다음 식에 의해 주어집니다
- (2) sj = \frac{-Cj(\mathbf{x}i)}{\nabla Cj \mathbf{M}-1 \nabla CjT} \quad
위치 변화량을 단순히 $k$로 스케일링함으로써 발생하는 바람직하지 않은 부작용은, 실질적인 유효 제약 강성(effective constraint stiffness)이 시간 단계(time step)와 수행된 제약 투영 횟수(number of constraint projections) 모두에 의존하게 된다는 점입니다. Müller 등 [2007]은 강성 계수의 지수적 스케일링(exponential scaling)을 통해 이 문제를 해결하려 시도했습니다. 하지만 이 방식은 시간 단계를 고려하지 않으며, 여러 제약 조건이 존재하는 경우 잘 정의된 해로 수렴하지 않습니다.
다음 섹션에서 우리는 잘 정의된 에너지 전위(energy potentials)와 직접적으로 대응되는 정규화된 제약 조건(regularized constraints)을 사용하는 기법을 개발하고, 이를 시간 단계와 반복 횟수에 독립적인 방식으로 해결하는 방법을 보여줄 것입니다.
Algorithm 1 XPBD simulation loop
1: predict position $\tilde{x} \Leftarrow x^n + \Delta t v^n + \Delta t^2 M^{-1} f_{ext}(x^n)$ 2: 3: initialize solve $x_0 \Leftarrow \tilde{x}$ 4: initialize multipliers $\lambda_0 \Leftarrow 0$ 5: while $i < solverIterations$ do 6: ㅤfor all constraints do 7: ㅤㅤcompute $\Delta \lambda$ using Eq (18) 8: ㅤㅤcompute $\Delta x$ using Eq (17) 9: ㅤㅤupdate $\lambda_{i+1} \Leftarrow \lambda_i + \Delta \lambda$ 10: ㅤㅤupdate $x_{i+1} \Leftarrow x_i + \Delta x$ 11: ㅤend for 12: ㅤ$i \Leftarrow i + 1$ 13: end while 14: 15: update positions $x^{n+1} \Leftarrow x_i$ 16: update velocities $v^{n+1} \Leftarrow \frac{1}{\Delta t}(x^{n+1} - x^n)$
4. Our Method
우리는 에너지 전위(energy potential) $U(\mathbf{x})$로부터 유도된 힘의 영향을 받는 뉴턴의 운동 방정식(Newton’s equations of motion)에서 시작하여 확장된 위치 기반 역학(XPBD) 알고리즘을 도출합니다.
- (3) \mathbf{M}\ddot{\mathbf{x}} = -\nabla UT(\mathbf{x}) \quad
여기서 $\mathbf{x} = [\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_n]^T$는 시스템 상태(system state)입니다. PBD에서 이는 흔히 단순한 입자 위치를 의미하지만, 강체 변환(rigid body transforms)과 같은 임의의 일반화된 좌표 모델(generalized coordinate model)을 나타낼 수도 있습니다. 기울기 연산자(gradient operator) $\nabla$는 편미분들의 행 벡터(row vector)라는 관례를 사용함에 유의하십시오.
우리는 운동 방정식 (3)에 대해 위치 수준의 암시적 시간 이산화(implicit time discretization)를 수행하며, 여기서 상첨자 $n$은 시간 단계 인덱스를 나타냅니다.
- (4) \mathbf{M} \left( \frac{\mathbf{x}n+1 - 2\mathbf{x}n + \mathbf{x}n-1}{\Delta t2} \right) = -\nabla UT(\mathbf{x}n+1) \quad
에너지 전위 $U(\mathbf{x})$는 제약 함수(constraint functions)의 벡터인 $\mathbf{C} = [C_1(\mathbf{x}), C_2(\mathbf{x}), \cdots, C_m(\mathbf{x})]^T$를 사용하여 다음과 같이 더 구체화될 수 있습니다.
- (5) U(\mathbf{x}) = \frac{1}{2} \mathbf{C}(\mathbf{x})T \boldsymbol{\alpha}-1 \mathbf{C}(\mathbf{x}) \quad
여기서 $\boldsymbol{\alpha}$는 강성의 역수(inverse stiffness)에 해당하는 블록 대각 컴플라이언스 행렬(block diagonal compliance matrix)입니다. 탄성 전위로부터 발생하는 힘은 $\mathbf{x}$에 대한 $U$의 음의 기울기(negative gradient)로 주어집니다.
- (6) \mathbf{f}elastic = -\nabla\mathbf{x} UT = -\nabla \mathbf{C}T \boldsymbol{\alpha}-1 \mathbf{C} \quad
우리는 Servin 등 [2006]의 방식을 따라, 라그랑주 승수(Lagrange multiplier)를 도입하여 힘을 방향과 스칼라 성분으로 분해함으로써 이를 유연한 제약 조건 공식(compliant constraint formulation)으로 변환합니다.
- (7) \boldsymbol{\lambda}elastic = -\tilde{\boldsymbol{\alpha}}-1 \mathbf{C}(\mathbf{x}) \quad
여기서 $\boldsymbol{\lambda}_{elastic} = [\lambda_1, \lambda_2, \cdots \lambda_m]^T$는 제약 승수(constraint multipliers)의 벡터입니다. 이후로는 명확성을 위해 필요한 경우가 아니면 하첨자를 생략하겠습니다. 식 (4)의 좌변에 있던 시간 단계를 컴플라이언스 행렬에 통합하여 $\tilde{\boldsymbol{\alpha}} = \frac{\boldsymbol{\alpha}}{\Delta t^2}$로 정의했음에 유의하십시오.
$\lambda$에 대한 식을 대입하면 다음과 같은 이산화된 제약 동역학 방정식을 얻습니다.
- (8) \mathbf{M}(\mathbf{x}n+1 - \tilde{\mathbf{x}}) - \nabla \mathbf{C}(\mathbf{x}n+1)T \boldsymbol{\lambda}n+1 = 0 \quad
- (9) \mathbf{C}(\mathbf{x}n+1) + \tilde{\boldsymbol{\alpha}}\boldsymbol{\lambda}n+1 = 0 \quad
여기서 $\tilde{\mathbf{x}} = 2\mathbf{x}^n - \mathbf{x}^{n-1} = \mathbf{x}^n + \Delta t \mathbf{v}^n$는 예측된 위치(predicted position) 또는 관성 위치(inertial position)라고 불립니다. 이 비선형 시스템을 해결하기 위해 우리는 뉴턴법에 기반한 고정점 반복(fixed point iteration)을 설계합니다. 이후 내용에서는 각 반복(iteration)별 미지수임을 강조하기 위해 하첨자 $(i + 1)$을 사용하며, 시간 단계 상첨자 $(n + 1)$은 생략합니다.
우리는 식 (8)과 (9)를 각각 $g$와 $h$로 정의합니다. 우리의 목표는 다음을 만족하는 $\mathbf{x}$와 $\lambda$를 찾는 것입니다.
- (10) g(\mathbf{x}, \boldsymbol{\lambda}) = 0 \quad
- (11) h(\mathbf{x}, \boldsymbol{\lambda}) = 0 \quad
식 (10)과 (11)을 선형화하면 다음과 같은 선형 뉴턴 서브문제(linear Newton subproblem)를 얻게 됩니다.
- (12) \begin{bmatrix} \mathbf{K} & -\nabla \mathbf{C}T(\mathbf{x}i) \\ \nabla \mathbf{C}(\mathbf{x}i) & \tilde{\boldsymbol{\alpha}} \end{bmatrix} \begin{bmatrix} \Delta \mathbf{x} \\ \Delta \boldsymbol{\lambda} \end{bmatrix} = - \begin{bmatrix} g(\mathbf{x}i, \boldsymbol{\lambda}i) \\ h(\mathbf{x}i, \boldsymbol{\lambda}i) \end{bmatrix} \quad
여기서 $\mathbf{K} = \frac{\partial g}{\partial \mathbf{x}}$입니다. 이 시스템을 통해 $\Delta \mathbf{x}$와 $\Delta \boldsymbol{\lambda}$를 구할 수 있으며, 이에 따라 위치와 승수(multipliers)가 업데이트됩니다.
- (13) \boldsymbol{\lambda}i+1 = \boldsymbol{\lambda}i + \Delta \boldsymbol{\lambda} \quad
- (14) \mathbf{x}i+1 = \mathbf{x}i + \Delta \mathbf{x} \quad
이는 $|\mathbf{x}_{i+1} - \mathbf{x}_i| \to 0$ 및 $|\boldsymbol{\lambda}_{i+1} - \boldsymbol{\lambda}_i| \to 0$을 만족하는 모든 수열에 대해 우리의 암시적 운동 방정식 (8, 9)를 만족시키는 고정점 반복(fixed-point iteration)입니다. 일반적으로 이 방법은 잘 작동하지만, 견고함(robustness)을 위해 라인 서치(line search) 전략이 필요할 수 있으며, 시스템 행렬을 계산하는 비용이 많이 들 수 있습니다. 특히 우리는 제약 조건의 헤시안(Hessians, 2차 도함수) 계산이 필요한 $\mathbf{K}$를 계산하는 것을 피하고자 합니다.
이제 우리는 구현을 단순화하고 기존 PBD와의 연결 고리를 만들어주는 두 가지 근사법(approximations)을 도입합니다. 첫째로, $\mathbf{K} \approx \mathbf{M}$이라는 근사를 사용합니다. 이는 기하학적 강성(geometric stiffness)과 제약 헤시안 항들을 생략하며, $O(\Delta t^2)$ 차수의 국소 오차를 발생시킵니다. 우리는 이 근사가 수렴 속도는 변화시킬 수 있지만, 고정점 반복의 전역 오차나 최종 해를 바꾸지는 않는다는 점에 주목합니다. 구체적으로, 이는 준 뉴턴법(quasi-Newton method)으로 볼 수 있습니다.
다음으로, 우리는 $g(\mathbf{x}_i, \boldsymbol{\lambda}_i) = 0$이라고 가정합니다. 이 가정은 $\mathbf{x}_0 = \tilde{\mathbf{x}}$ 및 $\boldsymbol{\lambda}_0 = 0$으로 초기화되는 첫 번째 뉴턴 반복에서 자명하게 참이라는 점을 통해 정당화됩니다. 또한, 제약 조건의 기울기(gradients)가 천천히 변한다면 이 값은 작게 유지될 것이며, 기울기가 일정할 때는 0이 될 것입니다. 나아가, 수정된 선형 시스템은 이제 현재 반복 지점인 $\mathbf{x}_i$에서 시작하여 제약 매니폴드(constraint manifold)로 향하는 질량 가중 최소화(mass-weighted minimization)의 최적 조건에 대응하게 됩니다.
이러한 근사들을 포함하여 업데이트된 선형 서브문제는 다음과 같이 주어집니다.
- (15) \begin{bmatrix} \mathbf{M} & -\nabla \mathbf{C}T(\mathbf{x}i) \\ \nabla \mathbf{C}(\mathbf{x}i) & \tilde{\boldsymbol{\alpha}} \end{bmatrix} \begin{bmatrix} \Delta \mathbf{x} \\ \Delta \boldsymbol{\lambda} \end{bmatrix} = - \begin{bmatrix} 0 \\ h(\mathbf{x}i, \boldsymbol{\lambda}i) \end{bmatrix} \quad
그다음 $\mathbf{M}$에 대해 슈어 보(Schur complement)를 취함으로써 미지수 $\Delta \boldsymbol{\lambda}$에 대한 다음과 같은 축소된 시스템을 얻을 수 있습니다.
- (16) \left[ \nabla \mathbf{C}(\mathbf{x}i) \mathbf{M}-1 \nabla \mathbf{C}(\mathbf{x}i)T + \tilde{\boldsymbol{\alpha}} \right] \Delta \boldsymbol{\lambda} = -\mathbf{C}(\mathbf{x}i) - \tilde{\boldsymbol{\alpha}} \boldsymbol{\lambda}i \quad
위치 업데이트는 다음 식을 계산함으로써 직접적으로 주어집니다.
- (17) \Delta \mathbf{x} = \mathbf{M}-1 \nabla \mathbf{C}(\mathbf{x}i)T \Delta \boldsymbol{\lambda} \quad
비록 본 방법이 제시하는 해가 암시적 운동 방정식을 정확하게 해결한다고는 할 수 없으나, 실제 적용 시 오차는 매우 작습니다. 우리는 섹션 6에서 본 방법의 정확성을 더 자세히 조사합니다.
💡 XPBD, 한마디로 무엇이 다른가?
기존 PBD가 "대충 이 정도 위치면 되겠지?"라고 눈대중으로 물체를 옮겼다면, XPBD는 "이 물체는 물리적으로 이만큼의 에너지를 가졌으니, 정확히 이만큼 움직여야 해"라고 수학적 근거를 가지고 계산하는 방식이다.
제시된 XPBD(Extended Position-Based Dynamics)의 핵심은 기존 PBD의 강성 문제를 해결하기 위해 에너지 전위($U$)와 라그랑주 승수($\lambda$)를 도입하여 물리적 엄밀함을 갖춘 것이다. 제약 조건을 탄성 에너지로 정의하고 이를 암시적 시간 이산화(Implicit Time Discretization)와 뉴턴법을 통해 풀어냄으로써, 시간 단계($\Delta t$)나 반복 횟수와 상관없이 일정한 강성을 유지하는 수치적 안정성을 확보했다.구현 측면에서는 복잡한 2차 도함수(Hessian) 계산을 피하기 위해 준 뉴턴법(Quasi-Newton) 근사를 사용하고, 대규모 행렬 식을 슈어 보(Schur Complement) 기법으로 축소했습니다. 이를 통해 최종적으로 각 제약 조건별로 위치와 승수를 독립적으로 업데이트하는 효율적인 구조를 도출했으며, 결과적으로 PBD의 계산적 이점(속도, 간결함)을 유지하면서도 정확한 물리적 구성 모델(Constitutive models)을 시뮬레이션할 수 있게 되었습니다.
4-1. A Gauss-Seidel Update
우리는 이제 선형 방정식 시스템 (16)에 대한 가우스-자이델(Gauss-Seidel) 해법을 고려함으로써 기존 PBD와의 연결 고리를 만듭니다. 인덱스 $j$를 가진 단일 제약 방정식을 취하면, 다음 식을 계산하여 해당 제약 조건의 라그랑주 승수 변화량($\Delta \lambda$)을 직접 구할 수 있습니다.
- (18) \Delta \lambdaj = \frac{-Cj(\mathbf{x}i) - \tilde{\alpha}j \lambdai,j}{\nabla Cj \mathbf{M}-1 \nabla CjT + \tilde{\alpha}j} \quad
이 식은 우리 방법의 핵심을 이룹니다. 제약 조건을 해결하는 동안, 우리는 먼저 단일 제약 조건에 대한 $\Delta \lambda_j$를 계산한 다음, 식 (13, 14)를 사용하여 시스템의 위치와 승수(multipliers)를 업데이트합니다. 알고리즘 1은 우리의 방법을 요약하고 있습니다. 이는 4행, 7행, 9행이 추가된 것을 제외하면 기존의 PBD 알고리즘과 동일합니다.
$\Delta \lambda_j$에 대한 식을 자세히 살펴보면, $\alpha_j = 0$인 경우 기존 PBD 알고리즘의 스케일링 인자 $s_j$ (식 2)와 정확히 일치함을 알 수 있습니다. 우리의 새로운 정의를 통해, 우리는 $s_j$를 무한히 단단한(infinitely stiff) 제약 조건에 대한 점진적인 라그랑주 승수의 변화량으로 이해할 수 있습니다.
컴플라이언스가 존재하는 경우($\alpha_j \ge 0$), 분자와 분모 모두에 추가적인 항이 나타납니다. 이 항들은 제약력을 제한하는 방식으로 제약 조건을 정규화(regularize)하며, 이는 식 (5)에서 주어진 탄성 전위(elastic potential)에 대응합니다.
식 (18)의 분자에는 $\lambda_{i,j}$가 포함되어 있습니다. 이는 현재 반복 $i$에서 제약 조건 $j$에 대한 총 라그랑주 승수(total Lagrange multiplier)입니다. 우리는 시스템 위치 외에도 이 변수를 저장하고 업데이트해야 합니다. 제약 조건당 하나의 스칼라 값을 추가로 저장하는 것은 적은 비용(overhead)이 들지만, 전체 제약력에 대한 유용한 정보를 제공합니다. 이는 파괴 가능한 조인트(breakable joints)와 같이 힘에 의존하는 효과를 구현하거나 햅틱 장치를 구동하는 데 사용될 수 있습니다.
5. Damping
우리의 방법은 역학적 에너지의 소산(dissipation)을 자연스럽게 포함하는 암시적 시간 단계 설정(implicit time stepping scheme)으로부터 유도되었습니다. 그럼에도 불구하고, 추가적인 제약 감쇠(constraint damping)를 모델링하는 것이 유용할 때가 있습니다. 이를 위해, 우리는 다음과 같이 레일리 소산 전위(Rayleigh dissipation potential)를 정의합니다.
- (19) D(\mathbf{x}, \mathbf{v}) = \frac{1}{2} \dot{\mathbf{C}}(\mathbf{x})T \boldsymbol{\beta} \dot{\mathbf{C}}(\mathbf{x}) \quad
- (20) = \frac{1}{2} \mathbf{v}T \nabla \mathbf{C}T \boldsymbol{\beta} \nabla \mathbf{C} \mathbf{v} \quad
여기서 $\boldsymbol{\beta}$는 제약 감쇠 계수(constraint damping coefficients)에 해당하는 블록 대각 행렬입니다. $\boldsymbol{\beta}$는 컴플라이언스($\alpha$)처럼 역수(inverse) 매개변수가 아니며, 감쇠 강성(damping stiffness)으로 설정되어야 한다는 점에 유의하십시오.
라그랑주 역학에 따르면, 소산 전위(dissipation potential)로부터 발생하는 힘은 속도($\mathbf{v}$)에 대한 $D$의 음의 기울기(negative gradient)로부터 유도됩니다.
- (21) \mathbf{f}damp = -\nabla\mathbf{v} DT = -\nabla \mathbf{C}T \boldsymbol{\beta} \nabla \mathbf{C} \mathbf{v} \quad
탄성(elastic)의 경우와 마찬가지로, 이 힘 벡터의 스칼라 성분을 제약 승수(constraint multiplier)로 분리합니다.
- (22) \lambdadamp = -\tilde{\beta} \dot{\mathbf{C}}(\mathbf{x}) = -\tilde{\beta} \nabla \mathbf{C} \mathbf{v} \quad
여기서 우리는 $\tilde{\beta} = \Delta t^2 \boldsymbol{\beta}$로 정의함으로써 감쇠 매개변수에 시간 단계를 포함시켰습니다. $\lambda_{damp}$를 단독으로 풀어서 감쇠 성분과 탄성 성분에 대한 제약력 승수를 각각 따로 구하는 것도 가능합니다. 하지만 대부분의 경우 우리는 전체 제약력만 알면 되며, 감쇠력은 탄성력과 동일한 방향으로 작용해야 하므로, 두 승수를 다음과 같이 하나의 방정식으로 결합할 수 있습니다.
- (23) \lambda = \lambdaelastic + \lambdadamp = -\tilde{\alpha}-1 \mathbf{C}(\mathbf{x}) - \tilde{\beta} \nabla \mathbf{C} \mathbf{v} \quad
식 (23)을 우리의 제약 조건 형식으로 재정렬하면 다음과 같습니다.
- (24) h(\mathbf{x}, \boldsymbol{\lambda}) = \mathbf{C}(\mathbf{x}) + \tilde{\alpha} \boldsymbol{\lambda} + \tilde{\alpha} \tilde{\beta} \nabla \mathbf{C} \mathbf{v} = 0 \quad
속도의 이산화 근사식인 $\mathbf{v} = \frac{1}{\Delta t}(\mathbf{x}^{n+1} - \mathbf{x}^n)$을 대입하고 $\lambda$에 대해 선형화하면, 다음과 같이 업데이트된 뉴턴 단계 식을 얻습니다.
- (25) \left[ \left(1 + \frac{\tilde{\alpha}\tilde{\beta}}{\Delta t}\right) \nabla \mathbf{C}(\mathbf{x}i) \mathbf{M}-1 \nabla \mathbf{C}(\mathbf{x}i)T + \tilde{\alpha} \right] \Delta \boldsymbol{\lambda} = -h(\mathbf{x}i, \boldsymbol{\lambda}i) \quad
단일 제약 방정식의 관점에서, 가우스-자이델(Gauss-Seidel) 업데이트 식은 이제 다음과 같이 주어집니다.
- (26) \Delta \lambdaj = \frac{-Cj(\mathbf{x}i) - \tilde{\alpha}j \lambdai,j - \gammaj \nabla Cj (\mathbf{x}i - \mathbf{x}n)}{(1 + \gammaj) \nabla Cj \mathbf{M}-1 \nabla CjT + \tilde{\alpha}j} \quad
$\tilde{\alpha}$와 $\tilde{\beta}$가 단순 대각 행렬인 일반적인 경우, $\gamma_j = \frac{\tilde{\alpha}_j \tilde{\beta}_j}{\Delta t}$가 됩니다. 즉, 이는 해당 제약 조건의 컴플라이언스와 감쇠 매개변수의 곱을 시간 단계로 스케일링한 값입니다. 우리는 여기서 추가된 모든 항이 쉽고 효율적으로 계산 가능하다는 점을 강조합니다.
위치에 저항하는 탄성력($\lambda_{elastic}$)과 속도에 저항하는 감쇠력($\lambda_{damp}$)을 하나의 라그랑주 승수($\lambda$)로 통합되었다.이를 통해 각 제약 조건 $j$에 대해 계산되는 라그랑주 승수 변화량 $\Delta \lambda_j$는 최종적으로 식 (26)에 의해 결정된다.
💡 그래서 PBD에서 어떻게 변경이 되었을까?1. 데이터 구조 추가 및 초기화- 라그랑주 승수 저장: 각 제약 조건(Constraint) 객체에 스칼라 변수 lambda_total을 추가한다.- 프레임별 초기화: 매 타임스텝 시작 시 lambda_total = 0으로 초기화하여 누적 제약력을 리셋한다.2. 제약 해항(Solver) 계산식 변경기존 PBD의 강성 계수($k$) 방식 대신 컴플라이언스($\alpha$) 기반 공식으로 교체한다.[변화량($\Delta \lambda$) 도출 공식]- \Delta \lambda = \frac{-C(\mathbf{x}) - \tilde{\alpha} \lambdatotal}{\sum w |\nabla C|2 + \tilde{\alpha}}
- 분자: 현재 위반량($-C$)에 누적 제약력($-\tilde{\alpha} \lambda_{total}$)을 반영하여 물리적 일관성을 유지한다.
- 분모: 컴플라이언스($\tilde{\alpha} = \alpha / \Delta t^2$)를 더해 수치적 안정성(정규화)을 확보한다.
3. 상태 업데이트 및 누적제약
충족 루프(Iteration) 내에서 위치와 힘을 동시에 갱신한다.
- 승수 누적: 계산된 변화량을 기존 승수에 즉시 더한다. (lambda_total += delta_lambda)
- 위치 보정: 별도의 강성 곱셈 없이 즉시 반영한다. ($\mathbf{x} \leftarrow \mathbf{x} + w \nabla C \Delta \lambda$)