FESADev/docs/hht-implicit-dynamics/formulation-inertia.md

# 유한회전 관성·초기화·HHT 상세

[주 정식화](formulation.md)의 구현 수식 상세다. 상태와 미결정 항목은 주 문서를 따른다. 아래 전개는 명시한 FESA 운동학에서의 자체 유도이며 Abaqus 내부 절차의 기술이 아니다.

2026-10-11 승인 [HHT-ADR-012](ADR.md)에 따라 이번 HHT 질량은 nodal lumped `T_L/M_L`다. Consistent 식은 비교와 미분 유도 이력으로 보존한다. 다른 physics·알고리즘·공개 output의 승인은 이 선택과 별개다.

## 1. 회전과 관성의 미분

`P=[phi]_×`, `theta=|phi|`에 대해

\[
J_l=I+\frac{1-\cos\theta}{\theta^2}P
+\frac{\theta-\sin\theta}{\theta^3}P^2,\qquad B=Dd=-[d]_\times J_l.
\]

0에서 scalar coefficient는 analytic series를 사용한다. `v_phi,a_phi`는 fixed-chart 시간미분이다. 방향 `p,s,r`에 대한 director 미분을
`B p=Dd[p]`, `C[p,s]=D²d[p,s]`, `D3[p,s,r]=D³d[p,s,r]`라 정의하면

\[
C[p,s]=(DJ_l[s]p)\times d+(J_lp)\times((J_ls)\times d),
\]
\[
\begin{aligned}
D3[p,s,r]={}&(D^2J_l[s,r]p)\times d+(DJ_l[s]p)\times B r\\
&+(DJ_l[r]p)\times B s+(J_lp)\times C[s,r].
\end{aligned}
\]

방향 교환에 대한 대칭성이 필요하다. 정확한 미분을 scalar coefficient의 비특이 series까지 구현하기 위한 동등한 정의는

\[
D^k\exp_P[H_1,\ldots,H_k]
=\sum_{m=k}^{\infty}\frac1{m!}
\sum_{\sigma\in S_k}\ \sum_{r_0+\cdots+r_k=m-k}
P^{r_0}H_{\sigma(1)}P^{r_1}\cdots H_{\sigma(k)}P^{r_k},\quad k=1,2,3.
\]

여기서 `H=[direction]_×`이고 오른쪽에 `d0`를 곱하면 필요한 director derivative다. 이는 수렴하는 해석적 정의이며 production에서 무한합을 그대로 계산하라는 뜻이 아니다. Rodrigues/Fréchet/정확한 미분으로 같은 도함수를 계산하고 finite difference는 독립 검증에 사용한다.

물질점에 대해 `H=D_q x`, `B_x[p,s]=D²_q x[p,s]`, `C_x[p,s,r]=D³_q x[p,s,r]`로 쓰면

\[
H p=\sum_iN_i(p_{ui}+zB_i p_{\phi i}),\quad
B_x[p,s]=\sum_iN_i z C_i[p_{\phi i},s_{\phi i}],
\]
\[
C_x[p,s,r]=\sum_iN_i zD3_i[p_{\phi i},s_{\phi i},r_{\phi i}],\quad
\dot x=Hv,\quad \ddot x=Ha+B_x[v,v].
\]

Consistent 관성 covector와 partial derivative action은 다음으로 직접 조립 가능하다:

\[
I_C=\int\rho H^T\{Ha+B_x[v,v]\}\,dV_0,
\quad M_C=\int\rho H^TH\,dV_0,
\]
\[
\begin{aligned}
D_qI_C[p]=\int\rho\big(&B_x[p,\cdot]^T\{Ha+B_x[v,v]\}\\
&+H^T\{B_x[p,a]+C_x[p,v,v]\}\big)dV_0,\\
D_vI_C[p]&=2\int\rho H^TB_x[v,p]\,dV_0,\qquad D_aI_C[p]=M_Cp.
\end{aligned}
\]

첫 항 `B_x[p,·]^T ẍ`를 빼면 단순 current-mass tangent에 불과하다. `B_x[p,a]`, `C_x[p,v,v]`, `2B_x[v,δv]`도 생략하지 않는다. `T`의 Euler–Lagrange 전개에서 `D H`의 대칭성에 의해 `d(H^T)/dt`와 `∂q T` 항이 상쇄되어 위 식을 얻는다. 원형상 밀도·measure는 시간에 따라 갱신하지 않는다.

채택한 nodal lumped 식에는 `x_i(z)=X_i+u_i+z d_i`, `H_i=D_q x_i`를 사용하여 위 적분을 `Σ_i ∫rho N_i (·)_i dV_0`로 바꾼다. `N_i`는 q와 무관하므로 동일한 derivative action이 성립한다. 별도의 경험적 gyroscopic 또는 Coriolis 보정항을 추가하지 않는다.

구체적인 현재 채택식은 다음과 같다. `B_i^x=D_q²x_i`, `C_i^x=D_q³x_i`이며
각 nodal mapping의 도함수는 해당 node의6개 성분에만 작용한다.

\[
T_L=\tfrac12\sum_i\int\rho N_i|H_iv|^2dV_0,\quad
M_L=\sum_i\int\rho N_iH_i^TH_i dV_0,
\quad I_L=\sum_i\int\rho N_iH_i^T\{H_ia+B_i^x[v,v]\}dV_0.
\]
\[
c_L=I_L(q,v,0),\quad
D_qI_L[p]=\sum_i\int\rho N_i\left((B_i^x[p,\cdot])^T\ddot x_i
+H_i^T\{B_i^x[p,a]+C_i^x[p,v,v]\}\right)dV_0,
\]
\[
D_vI_L[p]=2\sum_i\int\rho N_iH_i^TB_i^x[v,p]dV_0,
\qquad D_aI_L[p]=M_Lp.
\]

현재 위치·속도도 같은 nodal measure로 P/L을 계산한다. Geometry/Jacobian은 원래
bilinear shell 보간으로 검사하며 nodal kinetic mapping으로 geometry를 교체하지 않는다.
이를 통해 M만 바꾸고 기존 보간 기반 c/접선/에너지를 남기는 혼합 모델을 방지한다.

## 2. 채택 질량·비교식과 적분

주 문서의 scalar moment로

\[
(M_C)_{ij}=\begin{bmatrix}
C_{ij}^{(0)}I_3&C_{ij}^{(1)}B_j\\
C_{ij}^{(1)}B_i^T&C_{ij}^{(2)}B_i^TB_j
\end{bmatrix},\quad
(M_L)_{ii}=\begin{bmatrix}m_iI_3&s_iB_i\\s_iB_i^T&j_iB_i^TB_i\end{bmatrix},
\quad (M_L)_{ij}=0\ (i\ne j).
\]

`rho`의 차원은 `mass/length³`, `(C0,C1,C2)` 및 `(m,s,j)`의 차원은 `(mass,mass·length,mass·length²)`다. `M_uu/M_uphi/M_phiphi`는 각각 같은 세 차원, residual translation은 force, rotation은 force·length이다. `v_u,a_u`는 `length/time,length/time²`; `v_phi,a_phi,omega,omega_dot`는 `1/time,1/time²`. 각도 단위는 rad이며 SI를 추정하지 않는다.

Moments는 현행 [eigen2 mass 계약 M2–M3](../mitc4-shell-feast-eigen2/formulation-mass-stabilization.md#m2-consistent-scalar-moments와-적분)의 `3×3×3` Gauss–Legendre를 기준으로 하며, 그 scalar-moment lumping을 현재 HHT의 nodal kinetic measure에도 사용한다. 긴 유도는 선행 [mass 계약 M3](../mitc4-shell-feast-eigen/formulation.md#formulation-mass--m3-consistent-matrix와-충분한-quadrature)에 있다. 점은 `−sqrt(3/5),0,+sqrt(3/5)`, weight는 `5/9,8/9,5/9`, 순서는 증가하는 `xi→eta→zeta`다. Geometry가 bilinear degenerated shell이고 rho/t가 일정할 때 `Ni Nj z^k J0`의 축별 차수는 최대 4여서 이 moment 적분은 정확하다. `B_i`와 그 도함수는 nodal 값이므로 현재 관성과 tangent도 같은 moment를 사용한다.

모든 27 mass 점에서 finite positive `J0`와 기존 geometry predicate를 확인한다. Current trial도 내력 8점·tying·center·recovery와 추가 inertia 평가점에서 허용 geometry를 만족해야 한다. `abs(J)`, 점 생략, 강제 평면화는 없다. 물리 내력·물리 에너지와 section recovery는 기존 `2×2×2` 적분, drill potential은 기존 element-level 합을 보존한다.

Lumped positivity 조건은 `m_i>0`, `j_i−s_i²/m_i>0`다:

\[
2T_{L,i}=m_i\left|v_{ui}+(s_i/m_i)B_iv_{\phi i}\right|^2
+(j_i-s_i^2/m_i)|B_iv_{\phi i}|^2.
\]

두 질량식은 총질량 `Σm_i=∫rho dV0`, uniform translation energy, scalar thickness moments를 공유한다. 일반 속도장과 rigid angular motion의 orbital inertia는 다르다. 같은 현재 nodal `w_i=v_ui+z dDot_i`에 대해

\[
T_L-T_C=\frac14\int\rho\sum_{i,j}N_iN_j|w_i-w_j|^2dV_0\ge0.
\]

따라서 consistent와 lumped의 동적 응답이나 spectrum 동일성을 가정하지 않는다. 평판의 `s_i=0`, `j_i=t²m_i/12`를 warped geometry 일반식으로 사용하지 않는다. 사용자는 현재 nodal lumped kinetic mass를 선택했으며 기존 eigen의 lumped 계약은 변경하지 않는다.

## 3. 객관성·운동량·회복

공간에서 강체속도를 가하면 `v_ui=V+Omega×(Xi+ui)`, `v_phii=J_l^−1 Omega`이고 `dDot_i=Omega×d_i`다. Consistent는 `xDot=V+Omega×x`를 재현하고, lumped는 각 `x_i(z)`의 같은 강체속도를 적분한다. 따라서 각자의 운동에너지 measure로 회전 관성을 검증해야 한다.

Consistent의 물리적 총운동량·각운동량은

\[
P=\int\rho\dot x\,dV_0,\quad L_O=\int\rho(x-O)\times\dot x\,dV_0.
\]

Lumped는 `Σ_i ∫rho N_i xDot_i dV0`, `Σ_i ∫rho N_i(x_i−O)×xDot_i dV0`다. 이 정의에서 시간미분한 `Pdot,Ldot`와 외력·support reaction의 총 force/moment를 대조한다. Kinematic 변환은 algorithmic `A`가 아닌 물리적 `a`로 평가한다. Objective `U_phys+U_drill`의 전체 내력은 강체 가상운동에 일을 하지 않아야 한다. Numerical drill의 개별 nodal moment를 physical stress/resultant로 오인하지 않는다.

고정 proper global frame change 아래 모든 벡터·구속·하중을 함께 변환하면 scalar `T`, `U`와 virtual work가 불변이어야 한다. 기준축 일부만 회전시킨 비교는 이 검증에 해당하지 않는다. 중간 회전 branch를 통과하는 rechart는 지원하지 않는다.

## 4. Mass-null과 초기 일관성

`B_i`는 rank 2이며 instantaneous drilling 방향은

\[
n_i(q)=\begin{bmatrix}0\\J_l(\phi_i)^{-1}d_i\end{bmatrix},\quad
B_iJ_l^{-1}d_i=0.
\]

Consistent와 lumped 모두 이 방향에 `M n_i=0`, `n_i^T I=0`이다. 이 null은 일반적으로 `[0;d0]` 또는 고정 global `UR3`가 아니다. 필요한 geometry/rank 조건 아래 물리적 5 DOF/node 외 추가 null이 없는지 검사한다. 구속 후에는 `M_ff`와 free admissible variation 안에서의 null basis `N(q)`를 사용하며, constrained component를 단순 삭제한 nodal null을 그대로 basis라 부르지 않는다.

`g=f_int−F`라 두면 연속 대수 조건은

\[
\psi(q,t):=N(q)^Tg_f(q,t)=0.
\]

Physical energy는 director-parallel virtual work가 0이다. 본 force-only subset에서 `psi`는 numerical drill channel의 stationarity를 포함한다. Drilling stiffness는 algebraic coordinate를 결정할 수 있지만 mass rank를 바꾸지는 않는다. `N`의 q 미분도 포함하여

\[
\psi_q v_f+\psi_t=0,\qquad
\psi_q a_f+\psi_{qq}[v_f,v_f]+2\psi_{qt}v_f+\psi_{tt}=0
\]

를 초기 연속 일관성으로 요구한다. 첫 시간은 `0+`이며 이후 constant load의 우측 도함수 `Fdot=Fddot=0`를 사용한다. STEP jump를 impulse로 해석하여 초기 velocity를 바꾸지 않는다.

`q0=v0=0`, 고정 zero boundary에서는 다음 절차로 `a0`를 정할 수 있다.

1. `psi(q0,0+)=0`와 1차 조건을 확인한다. 여기서 물리적 initial acceleration을 풀기 전부터 불일관하면 실패다.
2. `b=F_f0−f_int,f(q0)−c_f(q0,v0)`에 `N0^T b=0`이 필요하다. 양의 질량 부분공간에서 `M_ff a_P=b`, `N0^T a_P=0`인 particular solution을 구한다.
3. `a_f0=a_P+N0 eta`를 놓고

\[
(\psi_qN0)\eta=-\psi_q a_P-\psi_{qq}[v0,v0]-2\psi_{qt}v0-\psi_{tt}
\]

를 푼다. 이 subset에서 우변 뒤 세 항은 0이다. `psi_qN0`가 nonsingular여야 초기 acceleration이 유일하게 정해진다. 단순 Moore–Penrose solve 결과를 그대로 가속도라고 채택하거나 massless acceleration을 0으로 놓지 않는다.

초기점에서는 `g0`가 병진 force를 포함하므로 `D N` 항을 형식적으로 삭제하지 않는다. 다만 zero rotation에서 실제 vanished 항은 식으로 확인할 수 있다. Null basis의 normalization/부호는 물리적 판정에 영향을 주지 않으며 source identity와 결정적 순서를 유지한다.

Bounds/regularity는 다음을 포함한다: 양의 `rho,t`; reference/current valid geometry; `||phi||<pi`와 기존 swing/twist 경계; 예상 mass rank의 국소 불변성; positive non-null inertia; nonsingular `psi_qN`; finite well-conditioned local branch. 시간 수렴 논의에는 해당 구간에서 이 positive/invertibility 조건이 균일하게 유지되어야 한다. 거리·rank·조건수 threshold의 수치는 Numerical Review 소유다. 실패를 small inertia, diagonal shift 또는 유령 drill torque로 보정하지 않는다.

시간적분에서는 endpoint `E_f=0`가 `psi=0`를 포함한다. 그러나 finite dt의 `psi_q v+psi_t` 및 2차 hidden-constraint defect가 자동으로 정확히 0이라는 주장은 하지 않는다. 이 quantity들을 평가하여 q/v/a의 시간 수렴과 함께 검증한다. Position residual이 작다는 사실만으로 massless velocity/acceleration의 안정성을 승인하지 않는다. 추가 projection/condensation이 필요하다는 결과가 나오면 시간 scheme 계약을 독립 검토하여 수정해야 한다.

## 5. 시간 정확도와 대안의 경계

### 5.1 일정 질량 HHT의 단순 확대를 배제하는 이유

`I=M(q)a+c(q,v)`, `g=f_int−F`, `b=c+g`라 하자. 잘못 확대하기 쉬운 식은

\[
M(q_{n+1})A_{n+1}+(1+\alpha)b_{n+1}-\alpha b_n=0,
\]

이며 Newmark에는 `A`를 사용한다. Smooth solution에서 필요한 보조변수의 전개는 `A=a+alpha h adot+O(h²)`다. 따라서 위 잔차의 leading correction은

\[
\alpha h(M\dot a+\dot b)=-\alpha h\dot M a
\]

로 일반적으로 0이 아니다. `c`를 넣고 exact tangent를 계산해도 이 오차는 없어지지 않는다. 이것은 nonlinear solver 수렴 오차가 아니라 시간 scheme의 문제다. 일정 질량 예제의 2차 결과는 이를 검출하지 못한다.

### 5.2 추천 후보의 근거와 자체 전개

[Brüls–Arnold, author postprint pp.3–6, 식 (5)–(10)](https://orbi.uliege.be/bitstream/2268/27088/1/pa_BrulsArnold2008_TheGeneralizedAlphaSchemeAsALinearMultistepIntegrator.pdf)는 endpoint 동적 평형과 별도 acceleration-like variable을 사용한다. 본 후보는 그 형식의 HHT parameter specialization이다. 해당 문헌의 일반론을 shell의 singular mass에 직접 대입하지 않는다.

**다음은 자체 Taylor 검산이다.** `A_n=(1+alpha)a_n−alpha a_(n−1)`이면 smooth 구간에서 `A_n=a_n+alpha h adot_n+O(h²)`다. Newmark velocity 식의 `h² adot` coefficient는

\[
\alpha+\gamma=1/2,
\]

로 정확하다. Position 식은 `h²a/2`를 재현하고 local defect는 `O(h³)`다. `A0=a0` startup의 velocity defect는 `O(h²)`로 시작값의 global 2차 요구에 맞는다. 보조 `A`를 없애면 fixed-step velocity recurrence는

\[
v_{n+1}-v_n=h\{-\alpha(1-\gamma)a_{n-1}
+[(1-\gamma)(1+\alpha)-\alpha\gamma]a_n
+\gamma(1+\alpha)a_{n+1}\}.
\]

세 coefficient 합은 1, `a_(n−1),a_n,a_(n+1)`의 시간 offset `−1,0,+1`에 대한 first moment는 `alpha+gamma=1/2`다. 여기에는 `Mdot` 잔여항이 없다. Endpoint equation이 nonsingular mass의 smooth `a=f(q,v,t)`를 정의하면 q/v의 2차 일관성과 안정성 검토를 할 수 있고 physical a도 같은 해의 정확도를 따른다. 이는 nonlinear mass의 smooth manufactured verification을 요구하는 근거이며 독립 수치 검토를 대신하지 않는다.

문헌의 nonsingular mechanical solve 및 regular index-3 조건과 달리 본 shell은 drilling direction에서 singular `M`을 가진다. §4의 `psi_qN` regularity와 constant rank의 실제 지원 범위, redundant full-coordinate 구현에서의 hidden-constraint drift·q/v/a 정확도는 **구현 후 필수 검증 항목**이다. 이 조건으로 reduced smooth dynamics를 설명할 수 있는 범위를 확인하되 모든 shell에 대한 보편 정리를 요구하지 않는다. [독립 검토](numerical-review.md#numerical-review)는 현 수식 후보를 대조했고, [독립 검산](numerical-review.md#numerical-review-probes)의 smooth variable-mass / moving-null 예제는 이를 지지한다. 이 제한된 evidence가 실제 shell/drilling 검증을 대신하지 않으며 유한회전 지원 범위를 축소하는 근거도 아니다.

### 5.3 다른 후보와 채택하지 않은 이유

| 후보 | 수식/장점 | 남은 충돌 |
| --- | --- | --- |
| Weighted mass | `[(1+alpha)M1−alpha M0] A1+(1+alpha)b1−alpha b0=0`; 위 `alpha h Mdot a` 항을 상쇄 | 서로 다른 current/old drilling null의 합은 null rank를 lift할 수 있다. Physical/algorithmic 의미와 DAE 정확도 확인 없이 채택하지 않음 |
| Ambient director | 독립 `y=(u_i,d_i)`에 일정 moment mass, `d_i·d_i=1` holonomic constraint와 drill orientation의 algebraic equation을 둔 HHT | Newmark ambient director rates의 constraint drift, rotation-vector state와의 정확한 rate mapping, additional multiplier/drill coordinate가 필요. 기존 6DOF state와 단순 동치라고 볼 수 없음 |
| Physical/algorithmic acceleration 분리 | 주 문서 후보; current `M_C(q)`만으로 physical endpoint equation을 풀며 null을 보존 | Shell mass-null의 시간 정확도·regularity 검증 필요. Consistent mass 선택이 이 시간적분 후보의 채택 승인까지 의미하지는 않음 |

기존 source는 nonlinear 내부력의 second derivative까지 계산하지만 관성 tangent에는 director의 third derivative가 필요하다. 기존 [small-rotation mass kernel](../../src/fesa/elements/mitc4_shell.cpp)의 moment는 재사용 가능한 기준이고 그 initial global mass를 유한회전 전체에 고정하는 것은 본 후보가 아니다. [정적 상태](../../include/fesa/analysis/nonlinear_analysis_state.h)는 `a/A` distinction을 제공하지 않으므로 현재 구현이 이 알고리즘을 지원한다고 보고하지 않는다.