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

# HHT implicit dynamics 정식화

- 상태: `independently-reviewed-candidate; physics-and-validation-pending`.
- Owner: `formulation-agent`, 2026-09-22. 문서 gate 산출물이며 구현 승인이나 수치 검증 통과가 아니다.
- 기준: [프로젝트 지침](../../AGENTS.md), [제품 PRD](../PRD.md), [공통 ADR](../ADR.md#adr-028), [Feature PRD](PRD.md), [요구사항](requirements.md), [연구](research.md#research).
- 관성의 전개식·초기 일관성·정확도 근거는 [formulation-inertia](formulation-inertia.md)가 소유한다. 영속 schema, C++ 저장 구조와 tolerance는 정의하지 않는다.

## 1. 결론과 적용 범위

**추천 후보는 유한회전 kinetic energy에서 유도한 관성력과, 물리적 가속도 `a` / 알고리즘 가속도 `A`를 분리하는 HHT 확장이다.** 일정 질량의 표준 HHT force weighting으로 정확히 환원한다. 단순히 `M(q)`를 기존 일정 질량 잔차에 대입하는 식은 추천하지 않는다. 이 차이는 [상세 §5](formulation-inertia.md#5-시간-정확도와-대안의-경계)에 유도한다.

**사용자가 shell 속도 보간의 운동에너지를 그대로 적분하는 consistent mass를 선택했다.** 따라서 HHT 기본값과 이번 구현 범위는 consistent mass다. 같은 운동학의 scalar-moment nodal lumped 식은 비교용 이론 및 기존 eigen 회귀 참고이며 이번 HHT 구현 요구가 아니다. 선행 eigen 결과나 동적 reference와의 적합도로 선택한 것이 아니다. 유한회전 지원 범위를 작은 회전으로 축소하지 않는다.

아래 HHT 확장과 `0+` 초기화는 FESA 채택 제안이다. 이번 consistent mass 선택이 다른 physics·알고리즘·공개 output 계약 전체의 승인을 의미하지 않는다.

| 구분 | 정식화 경계 |
| --- | --- |
| 물리 | 초기 형상 기준 Total Lagrangian FESA-MITC4, 큰 변위·유한회전, 작은 탄성 변형률, 균질 등방 선형 탄성, 고정 두께 |
| 시간 | 단일 step, `DIRECT`, `h=0.1`, `T=10`, `t_n=n h`, 초기 frame 및 100 increments. `INC=1000`은 increment 상한 |
| 하중/구속 | `t=0+`부터 node `2`의 global 1–3 force 각각 `−100000`; 외곽 지정 DOF 1–6은 모든 시점에 0 |
| 초기값 | `q_0=v_0=0`; `a_0`는 동적 초기 평형·대수 일관성으로 구하며 `A_0=a_0` |
| 제외 | 물리적 damping, contact, 소성·손상, follower/분포하중, nodal moment 입력 확대, 임의 amplitude, 시간변화 구속, adaptive/cutback, multi-step/restart, 다른 요소 dynamics |

원본 INP의 section-point 입력 `5`를 새 질량 또는 내력 quadrature 지정으로 해석하지 않는다. Abaqus 내부 요소·질량·회전·초기화와의 동등성은 목표가 아니다.

## 2. 좌표·운동학·재료와 내력

절점 좌표 순서는 정확히

\[
q_i=[u_x,u_y,u_z,\phi_X,\phi_Y,\phi_Z]_i^T,\quad v=\dot q,\quad a=\ddot q.
\]

`u`는 global translation, `phi`는 **초기 형상과 고정 global basis에 대한 total rotation vector**다. 각 절점에서

\[
R_i=\exp([\phi_i]_\times),\quad Q_i=R_iQ_i^0,\quad d_i=R_i d_i^0,
\quad \|\phi_i\|<\pi .
\]

Newton correction은 `q←q+Δq`이고, orientation은 그 total `phi`에서 다시 평가한다. Incremental rotation을 별도로 누적하거나 accepted step마다 chart/reference를 재설정하지 않는다. 기존 principal chart, director opposition, swing/twist, 양의 Jacobian 및 polar factor 유효성 경계를 보존한다. 실패 trial을 clamp하거나 재정규화하여 통과시키지 않는다.

\[
x(\xi,\eta,\zeta;q)=\sum_i N_i(X_i+u_i+z d_i),\quad z=t\zeta/2,
\quad dV_0=J_0\,d\xi d\eta d\zeta .
\]

`N_i`는 기존 bilinear MITC4 shape function이다. Interpolated director `ΣN_i d_i`를 단위화하지 않는다. 정확한 각속도와 각가속도는

\[
\omega_i=J_l(\phi_i)v_{\phi i},\qquad
\dot\omega_i=J_l(\phi_i)a_{\phi i}+DJ_l(\phi_i)[v_{\phi i}]v_{\phi i},
\]
\[
\dot d_i=\omega_i\times d_i,\qquad
\ddot d_i=\dot\omega_i\times d_i+\omega_i\times(\omega_i\times d_i).
\]

따라서 `v_phi`와 `a_phi`는 공간 각속도·각가속도가 아니다. `A_phi`는 이 물리 회복식에 넣지 않는다. `J_l`과 미분의 구현 가능한 정의는 [상세 §1](formulation-inertia.md#1-회전과-관성의-미분)에 있다.

물리 내력은 기존 [비선형 정식화 §§6–13](../nonlinear-static-newton-raphson/formulation.md)를 따른다. Green–Lagrange strain의 MITC assumed shear와 fixed-reference local engineering strain
`ē=[E11,E22,2E12,2E13,2E23]`에 대해

\[
U_{phys}=\frac12\int_{\Omega_0}\bar e^TC_5\bar e\,dV_0,\quad
C_5=\operatorname{diag}(C_{ps},(5/6)G I_2),\quad G=E/[2(1+\nu)],
\]
\[
C_{ps}=\frac{E}{1-\nu^2}
\begin{bmatrix}1&\nu&0\\\nu&1&0\\0&0&(1-\nu)/2\end{bmatrix}.
\]

Objective numerical drilling potential은 기존 polar-corotation / swing-twist angle `χ_ei`와 초기 고정 stiffness `k_de`를 그대로 사용한다:

\[
U_{drill}=\frac12\sum_{e,i}k_{de}\chi_{ei}(q)^2,\quad
f_{int}=\nabla_q(U_{phys}+U_{drill}),\quad K_{int}=\nabla_q^2(U_{phys}+U_{drill}).
\]

Drilling stiffness는 물리적 회전 관성이 아니다. `U_phys`와 `U_drill`을 구별하여 에너지 evidence를 계산하고, drill stress/resultant를 만들지 않는다. 재료 내부변수·소성 update는 없다. 재사용 대상 [evaluator](../../src/fesa/elements/nonlinear_mitc4_shell.cpp)의 `PhysicalEnergy`, `DrillingEnergy`, `EvaluateInternalResidualAndTangent`도 이 구분에 대조한다. 현재 공개 equilibrium/recovery 결과는 `physical_strain_energy`만 노출하므로 drill energy evidence의 전달은 후속 구현 작업이며 현 API가 이를 이미 제공한다고 보지 않는다.

## 3. 질량·관성·연속 평형

`B_i=D_phi d_i=−[d_i]_× J_l(phi_i)`이고 `C_ij^(k)=∫rho N_iN_j z^k dV_0`이다. 채택한 consistent `T_C`와 비교용 lumped `T_L`은 다음과 같다.

\[
T_C=\frac12\int_{\Omega_0}\rho\left|\sum_iN_i(v_{ui}+zB_iv_{\phi i})\right|^2dV_0,
\]
\[
T_L=\frac12\sum_i\left(m_i|v_{ui}|^2+2s_i v_{ui}\cdot B_iv_{\phi i}
+j_i|B_iv_{\phi i}|^2\right),\quad
(m_i,s_i,j_i)=\sum_j(C_{ij}^{(0)},C_{ij}^{(1)},C_{ij}^{(2)}).
\]

비교용 lumped 식은 scalar moment의 row sum 후 tensor를 재구성하며 global matrix row sum이 아니다. 이번 HHT에는 `T=T_C`, `M=M_C`를 사용하고 artificial drilling inertia는 허용하지 않는다. `T=½v^T M(q)v`의 질량행렬은 회전에 의존한다. 관성 covector는

\[
\boxed{I(q,v,a)=\frac{d}{dt}\partial_vT-\partial_qT=M(q)a+c(q,v)},
\qquad c_i=\sum_{j,k}\frac12(M_{ij,k}+M_{ik,j}-M_{jk,i})v_jv_k.
\]

고정 구속을 제외한 모든 admissible `δq`에 대해

\[
\delta q^T E=0,\quad E=I+f_{int}-F(t),
\qquad F(t)=F^\star\ (t\ge0+).
\]

이는 parent `rho_0 ẍ=Div_0 P`, `P=F_def S`의 가상일 이산화다. 강형의 `F_def`는 deformation gradient이며 nodal force `F`와 구별한다. 물리적 damping은 0이다. Tangent에는 `DM·a`, convective force의 configuration/velocity 미분을 모두 포함한다. 구체적인 matrix block, quadrature, derivative action과 null 조건은 [상세 §§1–4](formulation-inertia.md)에 고정한다.

## 4. 초기화와 HHT 시간적분 후보

음수 alpha 관례를 사용한다:

\[
-\frac13\le\alpha\le0,\quad \gamma=\frac12-\alpha,\quad
\beta=\frac14(1-\alpha)^2.
\]

기본값은 `alpha=−0.05`, `beta=0.275625`, `gamma=0.55`다. [공개 기본값 근거](research.md#research--3-reference-입력의-공개-의미)와 [Sandia의 결합계수](https://www.sandia.gov/files/sierra/SM_Users_5_20/main/implicit/integration.html)를 따른다. `alpha=0`은 평균가속도 Newmark다. 이 범위는 일반 유한회전 문제의 무조건 안정성이나 큰 시간간격의 정확도를 보증하지 않는다. 별도 `beta/gamma` 입력을 승인한 의미도 아니다.

**FESA 초기화 후보 계약:** frame 0은 `t=0+`, `q_0=v_0=0`, 전값 `F_0=F*`다. `M_ff(q_0)a_f0=F_f0−f_int,f(q_0)`와 [mass-null 일관성](formulation-inertia.md#4-mass-null과-초기-일관성)을 함께 풀어 `a_0`를 구한다. `A_0=a_0`, `g_0=f_int(q_0)−F_0`를 정의한다. `g_0=0` 또는 `a_0=0`을 임의 대입하지 않는다. 이는 공개 STEP 의미에 맞춘 FESA 계약 제안이며 Abaqus의 숨은 첫 history 생성 절차를 확인한 주장이 아니다.

[Brüls–Arnold (2008), 식 (5)–(8)](https://orbi.uliege.be/bitstream/2268/27088/1/pa_BrulsArnold2008_TheGeneralizedAlphaSchemeAsALinearMultistepIntegrator.pdf)의 보조가속도 형식을 `alpha_m=0`, `alpha_f=−alpha`로 전문화한다. 아래 `a`는 물리적 좌표 가속도, `A`는 시간적분 보조변수다:

\[
\boxed{A_{n+1}=(1+\alpha)a_{n+1}-\alpha a_n},
\]
\[
q_{n+1}=q_n+h v_n+h^2[(1/2-\beta)A_n+\beta A_{n+1}],
\qquad v_{n+1}=v_n+h[(1-\gamma)A_n+\gamma A_{n+1}],
\]
\[
\boxed{E_f(q_{n+1},v_{n+1},a_{n+1},t_{n+1})=0.}
\]

같은 fixed-global chart의 `a_n`만 history로 사용한다. 정확도 적용 조건과 source의 singular-mass 한계는 상세 §5에 명시한다. Shell drilling DAE까지 증명되었다고 보지 않는다.

Newton 미지수는 전체 좌표의 free 성분 `q_f,n+1`이다. `q_c=v_c=a_c=A_c=0`를 항상 유지한다. Predictor와 trial 종속량은

\[
q_p=q_n+h v_n+h^2(1/2-\beta)A_n,\quad v_p=v_n+h(1-\gamma)A_n,
\]
\[
A(q)=\frac{q-q_p}{\beta h^2},\quad
v(q)=v_p+\frac{\gamma}{\beta h}(q-q_p),\quad
a(q)=\frac{A(q)+\alpha a_n}{1+\alpha}.
\]

동일 trial에서 `E`를 평가하고 매 반복 full tangent를 갱신한다:

\[
\boxed{K_{eff}p=D_qI[p]+\frac{\gamma}{\beta h}D_vI[p]
+\frac{M(q)p}{(1+\alpha)\beta h^2}+K_{int}p},\quad
(K_{eff})_{ff}\Delta q_f=-E_f,\quad\Delta q_c=0.
\]

여기서 `D_qI`는 `v,a`를 고정한 편미분이다. 외력은 고정 global force이므로 `D_qF=0`이다. 관성 convective 항 때문에 `K_eff` 대칭을 일반적으로 가정하지 않는다. 모든 DOF가 구속되면 free solve는 비어 있고 full reaction을 회복한다. Singular mass만으로 모델을 거부하거나 epsilon mass를 추가하지 않으며, 유효 접선의 singularity·대수 regularity 실패는 별도로 진단한다.

## 5. 잔차·반력·에너지의 구별

일정 `M`에서 `g=f_int−F`라 두면, 앞선 endpoint physical equilibrium `M a_n+g_n=0`을 이용하여

\[
R_{HHT}:=M A_{n+1}+(1+\alpha)g_{n+1}-\alpha g_n
=(1+\alpha)E_{n+1}-\alpha E_n.
\]

따라서 free 방정식은 표준 force-weighted HHT와 같다. **표준 HHT의 Newmark 가속도를 이 문서에서는 `A`라 부른다.** 이를 `a`로 오인하면 출력·반력이 달라진다. 구성 의존 질량에는 위 `M A+weighted g` 식을 그대로 적용하지 않는다. 추천 후보의 Newton residual은 endpoint `E_f`이며, HHT filter 결함 `A_{n+1}−(1+alpha)a_{n+1}+alpha a_n`은 별도 운동학 evidence다.

구속 selector `C_u,C_phi`에 대한 endpoint 반력은

\[
F_i^{react}=C_{u,i}E_{u,i},\quad m_{\phi,i}^{react}=C_{\phi,i}E_{\phi,i},\quad
M_i^{react}=J_l(\phi_i)^{-T}m_{\phi,i}^{react}.
\]

부호는 structure에 작용하는 support force/moment다. Free `E_f`는 평형 evidence이고 반력이 아니다. `R_HHT,c`는 이전 반력이 섞인 알고리즘량이며 endpoint 반력으로 저장하면 안 된다. Moment는 `m_phi·δphi=M_spatial·δtheta`에 의해 변환한다. Numerical drill의 force/work는 물리적 section stress와 분리한다.

연속계의 에너지 관계는 `d(T+U_phys+U_drill)/dt=v^TF`다. 일정 dead force와 고정 경계에서는

\[
W_{0\to n}=F^{\star T}(q_n-q_0),\quad
\Delta\mathcal E_n=T(q_n,v_n)+U_{phys}(q_n)+U_{drill}(q_n)-\mathcal E_0-W_{0\to n}.
\]

`ΔE_n`은 이산 에너지 결함이며 물리 damping energy가 아니다. HHT의 일반 비선형 문제에서 매 frame 단조 감소나 정확한 momentum 보존을 요구하지 않는다. 무감쇠 선형 `alpha=0`의 보존성, 음수 alpha의 선형 고주파 소산과 비선형 work/energy 수렴은 따로 검증한다.

## 6. 증분·실패·회복 순서

1. 승인 subset, source identity, 재료·geometry·quadrature 및 고정 구속을 확인한다. Consistent mass의 moment·rank/positivity 조건과 초기 대수 일관성을 검사한다.
2. 초기 상태 및 `F_0`, `a_0`, `A_0`를 완성하여 frame 0을 검증한다. 초기화 실패를 zero acceleration으로 대체하지 않는다.
3. `t_n=n h`의 직전 committed `q,v,a,A`에서 predictor를 만들고, 고정 구속을 적용한다. `INC`와 정확한 종료 increment를 검사한다.
4. 각 Newton trial에서 total rotation, 내력·tangent, 관성·tangent와 `E`를 같은 state로 계산하고 free correction을 푼다. 병진/회전 correction 및 force/moment residual을 각각 scale-aware로 검사한다. 수치 한계·iteration 상한은 Numerical Review/I/O 소유다.
5. 수렴 후 같은 trial에서 physical acceleration, triad/director, endpoint reaction, 물리적 strain/stress/resultant, 에너지 및 null-constraint evidence를 회복·검증한다. 고정 reference frame과 element Gauss/section-point identity는 기존 정적 회복을 보존하며 nodal averaging은 하지 않는다.
6. `time,increment,q,v,a,A,orientation,forces,recovery,evidence`가 한 identity를 이루는 candidate만 원자적으로 commit한다. 모든 101 accepted frames를 보존한다. Newton iteration은 time frame이 아니다.
7. 비수렴·invalid chart/geometry·nonfinite·선형 solve·회복·출력 실패는 trial을 폐기하고 마지막 committed 상태/history와 기존 최종 HDF5를 보존한다. 고정 `DIRECT`를 cutback/adaptive로 바꾸거나 실패 frame을 성공으로 내보내지 않는다.

공식 출력은 lexical sibling `<stem>_results.h5`다. `U/UR`, 물리적 `v/a`, `omega/omega_dot`, algorithmic `A`, 반력과 에너지의 **공개 여부 및 영속 schema는 미결정**이다. Algorithmic `A`를 physical acceleration이라는 이름으로 공개해서는 안 된다. 출력 선택과 무관하게 검증에는 서로 구별되는 내부 quantity가 필요하다.

## 7. 검증 인계와 미해결 결정

| 항목 | 필수 evidence / 남은 경계 |
| --- | --- |
| HHT 계수·initialization | 표준 일정 질량 HHT와 algebraic reduction; SDOF 자유진동/STEP onset 및 `a_0=F/m`; `alpha=0` 한계 |
| 관성·회전 | Consistent mass의 총질량, uniform translation, 유한회전 angular motion, `d·ḋ=0`, `d·d̈+|ḋ|²=0`, objectivity, kinetic-energy derivative·tangent directional derivative |
| mass-null | 구성 의존 null, `ψ=0`과 초기 1·2차 일관성; coupled massless DAE에서 q/v/a 시간 수렴, finite-rotation drill coupling. 이 source theorem의 직접 적용은 미확인 |
| 정확도 | Smooth `M(q)` manufactured case의 q/v/a를 h,h/2,h/4로 검산; constant-M와 variable-M를 분리. STEP의 시작 불연속과 이후 smooth 구간을 분리 |
| 에너지·반력 | 관성 구속 평형, generalized/spatial moment work, algorithmic/endpoint reaction 구별, `T/U_phys/U_drill/W`와 momentum evidence |
| 구현 lifecycle | 동일 state tangent, 실패 rollback, time/frame 순서, deterministic assembly와 반복 결과 |
| 외부 reference | 지정 원본 101 frames × 9 nodes, frame/instance/node/component의 exact mapping; U blocking / UR warning-only 및 공통 TOLERANCE. V/A/energy correctness는 해당 CSV에 coverage 없음 |

수식 후보의 [독립 검토](numerical-review.md#numerical-review)는 완료되었으며 표시된 수식에서 새 오류는 발견되지 않았다. 질량 선택은 consistent로 확정되었다. 구현 전 남은 계약은 (1) HHT 확장·초기화의 FESA 채택과 나머지 physics 범위, (2) initial algebraic regularity와 trial 실패 판정의 수치 기준, (3) 공식 output 의미의 공개 범위/I/O다. 이 선택 후 구현할 수 있는 수식 후보이며, 실제 유한회전 shell의 drilling DAE·관성 tangent·q/v/a 시간 수렴은 **구현 후 기능 acceptance에 필요한 검증**이다. 부모 문헌의 nonsingular-mass 조건을 shell mass-null에 자동 적용하지 않되, 모든 shell DAE의 보편적 수렴 증명을 코드 작성의 선행 조건으로 요구하지 않는다. 현재 문서 완료나 제한된 probe를 전체 solver 승인·검증 통과로 표시하지 않는다.

이번 owner 작업은 기존 정식화·source 대조, 일차 문헌 대조와 수식 유도 및 문서 정적 검토다. [독립 시간정확도 검산](numerical-review.md#numerical-review-probes)은 variable-mass SDOF와 moving-null 예제에서 추천 후보의 약 2차 q/v/physical-a 수렴을 관찰했고 naive 확대의 약 1차 오차를 구별했다. 이는 shell solver 검증이나 일반 drilling DAE 증명을 대신하지 않는다. Build/test·solver/reference 실행은 문서 작업 범위가 아니므로 미실행이다. 독립 검산과 Numerical Review 결과는 해당 owner의 보고서가 소유하며 본 문서가 스스로 승인하지 않는다.