FESADev/docs/hex8-structural-analysis/formulation.md

# Hex8 정식화와 검증 불변량

> 현재 구현·검증은 [기능 요약](README.md)을 먼저 따른다. 아래 준비 상태·owner 인계와
> 옛 source 문서명은 작성 당시 이력이며, 원문은 Git 이력에서 확인한다.

적용: [PRD](README.md)의 FEM-REQ-HEX8-002–007,012–018,021.
물리 선택은 [ADR](README.md)에 있다. 이 문서는 승인된 식의 구현용 유도이며
외부 결과에 맞추어 physics를 바꾸는 권한이 아니다. 수치 controls는
[별도 원본](README.md), 출력 tensor/identity는 [I/O](README.md)를 따른다.
2026-10-09 승인된 Solid HHT 내부 정밀도 변경의 구현 계약은
[precision formulation](precision-formulation.md)이다. 동적 paired 경로의 상대변위 gradient,
전체 상태·계수 연산과 반올림 출력 의미는 해당 문서를 함께 적용한다.
사용자 선택 B로 승인된 현재 모델은 **8점 deviatoric 적분 + 기준 체적가중 평균 Green trace
energy**다. 선형/비선형 정적, 무응력 원형상 frequency, 선형/비선형 HHT 모두 같은 energy의
선형극한·내력·Hessian을 사용한다. 이전 비수정 pointwise full8의 수식/evidence는 변경 전
이력이며 현재 모델의 PASS로 승계하지 않는다. Fbar·평균 determinant·log-strain law는 아니다.

## Shape·기준 geometry·적분

절점 signs는 `---,+--,++-,-+-,--+,+-+,+++,-++`다. Point signs는
`---,+--,-+-,++-,--+,+-+,-++,+++`, coordinate는 각 ±1/√3이고 weight1이다.

\[
N_a={1\over8}(1+\xi_a\xi)(1+\eta_a\eta)(1+\zeta_a\zeta),\quad
J_{0,iJ}=\sum_aX_{ai}N_{a,J},\quad G_a=J_0^{-T}\nabla_\xi N_a.
\]

J0 column은 자연 ξ/η/ζ 미분이다. Gauss reference 위치는 ΣNaXa, volume weight는 det J0다.
Translation 차이를 먼저 사용해 large coordinate offset의 취소 오차를 줄일 수 있지만
기준 좌표와 source identity는 변경하지 않는다. Shape의 partition of unity·linear reproduction와
gradient 합0을 검증한다. Invalid det/inverse는 abs/clamp나 corner 재배열로 숨기지 않는다.

실제 상대변위 계산과 그 미분을 모든 절차에서 일치시킨다. Local 절점 index는0–7이다.
기존 immutable F64 cache의 \(G_{ag},w_g\)를 그대로 사용하고

\[
 H_g=\sum_{a=1}^7(u_a-u_0)\otimes G_{ag},\quad
 \widehat G_{ag}=G_{ag}\ (a>0),\quad
 \widehat G_{0g}=-\sum_{a=1}^7G_{ag},\quad
 D_g[p]=\sum_{a=0}^7p_a\otimes\widehat G_{ag}.
\]

\(D_g[p]\)는 방향 p에 대한 \(\delta F_g=\delta H_g\)다. F64 static/frequency 및 paired
dynamic은 각 경로의 산술로 같은 순서의 coefficient 합을 평가한다. Relative H를 사용하고
force/tangent에는 별도의 cached \(G_{0g}\)를 쓰는 혼합은 허용하지 않는다.
새 quadrature로 coefficient를 재생성하지 않으며 coefficient/합산 roundoff는 독립 검증한다.

## Cell 전체의 Jacobian 인증

Trilinear map의 det J는 각 자연 좌표에 대해 최대2차인 tensor polynomial이다.
`s=(ξ+1)/2,t=(η+1)/2,r=(ζ+1)/2`에서 degree(2,2,2) Bernstein basis로 표현할 수 있다.
각 축 -1/0/+1의27개 determinant 값에서 다음1D 변환을 세 축에 적용한다.

\[
b_0=p(-1),\quad b_2=p(1),\quad b_1=2p(0)-(p(-1)+p(1))/2.
\]

Bernstein basis는 nonnegative이고 합1이므로 모든 coefficient의 양의 lower bound는
cell 전체 det>0를 보증한다. 모든 sample가 양수인 것만으로는 충분하지 않다.
Coefficient bound가 불확실하면 de Casteljau midpoint subdivision 후 child cell을 검사한다.
Sample의 nonpositive upper bound 또는 완전히 nonpositive cell은 invalid이며, depth/resource
한계에서도 positive를 인증하지 못하면 uncertain geometry failure다.

인증용 연산은 outward-rounded interval arithmetic으로 수행한다. Addition/subtraction/
multiplication·변환·subdivision·current Xa+ua까지 실제 FP 오차를 포함하며 NaN/Inf/overflow를
거부한다. 단순 coefficient estimate에 임의 epsilon을 더해 positive라고 선언하지 않는다.
Current geometry에도 같은 인증을 적용하고 point det F>0를 확인한다. Trial failure가 reference
cache나 accepted state를 변경하지 않는다. 세부 종료 controls는 numeric-controls를 따른다.

## 재료·체적 평균 energy·선형 K0

\(\mu=E/[2(1+\nu)]\), \(\lambda=E\nu/[(1+\nu)(1-2\nu)]\),
\(\kappa=\lambda+2\mu/3=E/[3(1-2\nu)]\)는 실수 수학 identity다.
실행 coefficient는 \(\nu\to-1\)에서 \(\lambda+2\mu/3\)의 취소를 피하도록
**E/ν에서 직접 F64 \(\kappa=E/[3(1-2\nu)]\)를 한 번 계산**하고 finite positive를 검증한다.
계산 순서는 `2*nu`, `1-(2*nu)`, `3*(1-(2*nu))`, `E/denominator`의 RN64 연산이다.
`E/3` 선계산 등으로 재배열하지 않는다. Overflow/nonfinite 및 underflow-to-zero κ는
지원 산술 실패이며 값 보정은 없다. 기존 λ 조합의 bulk guard로 유효한 직접 κ를 다시 거부하지 않는다.
현재 실행 coefficient는 기존 정의의 F64 μ와 직접 F64 κ다. λ는 위 실수 identity와 변경 전
이력을 설명하는 기호이며 새 law의 별도 실행 cache/guard가 아니다. Paired 경로는 μ/κ를 exact lift한다.
변경 전 rounded μ/λ로 다시 만든 κ와의 bit equality를 요구하지 않으며 near-bound coefficient 오차를
독립 검토한다. \(E>0,-1<\nu<.5\)와 실제 positive finite shear/bulk를 요구하며 새로운 ν 범위
제한이나 clamp는 없다. λ 조합이나 `3*kappa` 등의 파생 overflow로 material 범위를 추가로
제한하지 않는다. 실제 element/trial 연산 graph가 nonfinite/지원 범위 실패이면 해당 평가를 실패시킨다.
\(\operatorname{dev}T=T-\operatorname{tr}(T)I/3\)다.
\(w_g=\det J_{0g}>0\), \(V=\sum_{g=1}^8w_g>0\)는 불변 기준 coefficient다.
Deviatoric 항·평균의 numerator/denominator 모두 동일 2×2×2 점을 사용한다.
평균은 기준 volume-weighted 값이며 center strain·단순8점 평균·current-volume 평균으로 바꾸지 않는다.

선형은 실제 \(\epsilon_g=\operatorname{sym}H_g\)에서

\[
 \bar\theta_{lin}=\sum_gw_g\operatorname{tr}\epsilon_g/V,\quad
 U_{lin}=\mu\sum_gw_g\operatorname{dev}\epsilon_g:\operatorname{dev}\epsilon_g
                        +\kappa V\bar\theta_{lin}^2/2,
\]
\[
 \sigma_g=2\mu\operatorname{dev}\epsilon_g+\kappa\bar\theta_{lin}I,
 \quad f_{ai}=\sum_gw_g\sigma_{g,iJ}\widehat G_{ag,J},\quad
 f=K_0u,\quad U_{lin}=u^TK_0u/2.
\]

Engineering strain vector는 `[ε11,ε22,ε33,2ε12,2ε13,2ε23]`이고 stress의 shear는 실제 값이다.
\(Q_g^{ai}=\operatorname{sym}(e_i\otimes\widehat G_{ag})\),
\(b^0_{ai}=\sum_gw_g\widehat G_{ag,i}\)이면 zero-state Hessian은

\[
 (K_0)_{ai,bk}=2\mu\sum_gw_g\operatorname{dev}Q_g^{ai}:\operatorname{dev}Q_g^{bk}
                             +\kappa b^0_{ai}b^0_{bk}/V.
\]

첫 항의 point block은 \(\mu\delta_{ik}\widehat G_a\cdot\widehat G_b+
\mu\widehat G_{a,k}\widehat G_{b,i}-(2\mu/3)\widehat G_{a,i}\widehat G_{b,k}\)다.
마지막 체적 외적은 모든 GP의 coupling을 포함한다. 기존 비수정 point-local
\(\sum B^TCB w\)를 선형 정적·frequency·linear HHT의 K0로 재사용하지 않는다.

## Total Lagrangian energy·내력·전체 Hessian

\[
 F_g=I+H_g,\quad E_g=\operatorname{sym}H_g+H_g^TH_g/2=(F_g^TF_g-I)/2,
 \quad\theta_g=\operatorname{tr}E_g,\quad\bar\theta=\sum_gw_g\theta_g/V,
\]
\[
 U=\mu\sum_gw_g\operatorname{dev}E_g:\operatorname{dev}E_g
                                      +\kappa V\bar\theta^2/2,
 \quad S_g=2\mu\operatorname{dev}E_g+\kappa\bar\theta I,\quad P_g=F_gS_g.
\]

이 식은 요소의 discrete nonlocal energy다. Homogeneous material과 reference coefficient 외에
소성/history state는 없다. 매 trial의8점 E를 평가한 뒤 같은 trial의 평균/S/P를 계산한다.
\(Z_g[p]=\operatorname{sym}(F_g^TD_g[p])\),
\(Y_g[p,q]=\operatorname{sym}(D_g[p]^TD_g[q])\)라 두면

\[
 \delta U[p]=\sum_gw_gS_g:Z_g[p],\quad
 f_{ai}=\sum_gw_gP_{g,iJ}\widehat G_{ag,J},
\]
\[
 \delta^2U[p,q]=\sum_gw_g\{2\mu\operatorname{dev}Z_g[p]:\operatorname{dev}Z_g[q]
                                      +S_g:Y_g[p,q]\}
 +\frac{\kappa}{V}\Big(\sum_gw_g\operatorname{tr}Z_g[p]\Big)
                  \Big(\sum_hw_h\operatorname{tr}Z_h[q]\Big).
\]

\(Z_g^{ai}=\operatorname{sym}(F_g^T(e_i\otimes\widehat G_{ag}))\),
\(b_{ai}=\sum_gw_g(F_g\widehat G_{ag})_i\)를 사용한 구현용 block은

\[
 (K_{mat})_{ai,bk}=2\mu\sum_gw_g\operatorname{dev}Z_g^{ai}:\operatorname{dev}Z_g^{bk}
                                        +\kappa b_{ai}b_{bk}/V,
\]
\[
 (K_{geo})_{ai,bk}=\delta_{ik}\sum_gw_g\widehat G_{ag}^TS_g\widehat G_{bg},
 \quad K_t=K_{mat}+K_{geo}.
\]

\(\delta S_g=2\mu\operatorname{dev}\delta E_g+\kappa\delta\bar\theta I\)이며
\(\delta\bar\theta=\sum_hw_h\operatorname{tr}\delta E_h/V\)다.
체적 Hessian은 \(g\ne h\)도 포함한다. 평균을 Newton 동안 고정하거나 GP별 기존 StVK C로
접선을 만들거나 stress만 평균화하면 이 energy의 내력·Hessian이 아니다.
Kt는 symmetric energy Hessian이나 stressed trial에서 SPD라고 가정하지 않는다.
Load-independent CLOAD tangent는0이고, zero state에서 F=I/S=0/Kgeo=0/Kt=K0다.

### 객관성·rank·복구 의미

실수 연산과 exact geometric completeness \(\sum_a\widehat G_{ag}=0\),
\(\sum_aX_a\otimes\widehat G_{ag}=I\)에서 \(u_a=(R-I)X_a+c\)이면 F=R/E=0/U=f=0이다.
Cached F64 coefficient/gradient 합/paired 반올림의 실제 rigid rotation은 coefficient와 연산
roundoff의 bounded 검사 대상이다. Gradient partition sum만으로 affine reproduction/exact0을
주장하지 않는다. 이미 주어진 F를 global R로 회전하면 E/평균/U는 invariant이고 force/tangent는
covariant다. 이 객관성은 평균 Green trace에 관한 것이며 determinant 평균화가 아니다.

K0는 nonnegative quadratic energy의 Hessian이다. 유효 brick의 local6 rigid modes/rank18,
warped cell의 추가 영모드·conditioning 및 실제 BC 뒤 mechanism/singularity를 검증한다.
특정 brick rank를 모든 admissible cell의 엄격 coercivity 증명으로 대신하지 않는다.
\(U_{old}-U=\kappa\sum_gw_g(\theta_g-\bar\theta)^2/2\ge0\)에서 Uold는 변경 전 비수정
pointwise full8의 **같은 μ/κ로 쓴 수학적 energy**다. 이전 실제 λ/μ rounded cache와 현재 직접 κ의
bit equality를 이 항등식에서 추론하지 않는다. 체적변형률 변동의 energy가 감소하므로
새 잠김·공간수렴·rank 검사가 필요하다.
Affine F에서는 위 같은 μ/κ/gradient로 정의한 Uold의 U/force/stress와 수학적으로 일치하지만
이전 rounded λ/raw G0 경로와의 byte equality나 모든 비affine direction의 Hessian 일치를 뜻하지 않는다.

공식 비선형 stress는 실제 \(\sigma_g=F_gS_gF_g^T/\det F_g\), strain은 **실제 기준 global Green E**다.
S의 isotropic 항이 요소 내 공통이어도 finite Cauchy pressure가 반드시 상수는 아니다.
선형 공식 strain은 실제 infinitesimal tensor이며 stress는 위 평균 trace law다.
8점 `[11,22,33,12,13,23]` 실제 tensor shear를 그대로 저장하고 평균/외삽/nodal smoothing하지 않는다.
저장 strain을 \(\operatorname{dev}E+\bar\theta I/3\)로 치환하지 않는다. Point stress는 같은 point의
strain만으로 결정되는 이전 pointwise law가 아니라 element mean에 의존한다.

## Reference row-sum lumped mass와 generalized eigen

\[
(M_C)_{ai,bj}=\delta_{ij}\sum_g\rho_0N_aN_bw_g,
\quad m_a=\sum_b(M_C)_{ai,bi}=\sum_g\rho_0N_aw_g,
\quad (M_L)_{ai,bj}=\delta_{ab}\delta_{ij}m_a,
\quad T={1\over2}v^TM_Lv.
\]

M은 reference-state 상수다. Positive weights와8점 interpolation matrix의 full rank에서 SPD가
따른다. Axis-aligned brick의 scalar block은 \(\rho V/216\)에 각 축 sign이 같으면2/다르면1인
세 factor의 곱이다. Direction별 총합은 rho V, row sum은 rho V/8다. 2026-10-11
[HEX8-ADR-010](README.md) 승인에 따라 production은 full24×24 M_C의 행 합 M_L을 사용한다.
행 합은 BC restriction 전에 계산한다. Warped 요소에서는 일반적으로 m_a가 같지 않다.
Positive rho/N/weights에서 m_a>0이며 M_L은 diagonal SPD다.

같은 nodal v에서 `T_L−T_C=0.25*sum_g rho*w_g*sum_ab N_a*N_b*|v_a−v_b|²>=0`이다.
따라서 M_L−M_C는 positive semidefinite이며 다음 consistent lower bound는 M_L에도
보수적으로 유효하다. Uniform translation energy와 총질량을 보존하지만 일반 속도장의
에너지·orbital inertia·고유치는 달라진다. Shell의 finite-rotation nodal block 방식과 구별한다.

8점 shape evaluation의 minimum singular value 제곱은1/27이다. Exact local quadrature에서
\(m_{e,lower}=\rho\min_gw_g/27\)는 local M의 positive lower bound다. Global scatter·principal free
restriction에도 active-node incidence를 고려한 lower bound가 보존된다. 예를 들어 모든 active
node에 적어도 하나의 element가 있는 경우 min_e bound는 conservative global lower bound다.
실제 FP matrix의 기존 roundoff guard에는 full row-sum 누적/조립 graph를 포함해 검토하고
SPD evidence를 확인해야 한다. 이 수학적 식만으로 설치 backend의
실행을 검증했다고 주장하지 않는다.

무응력 원형상에서 위 새 energy의 K0를 사용해 \(K_{ff}\phi=\lambda M_{ff}\phi\)를 풀며 constrained component0으로 복구하고
\(\phi^TM\phi=1\), \(f=\sqrt\lambda/(2\pi)\)를 구분한다. Proven mass lower bound가 있다면
\(\lambda_{max}\le\|K_{ff}\|_\infty/m_{lower}\)는 search upper-bound 근거가 된다.
Diagonal ratio classification scale 자체는 upper bound가 아니다. 검색/count의 실제 증거와
모든 backend pair의 original K/M 잔차·mass norm/Gram을 검증한 뒤 최저 N을 공개한다.
Repeated modes는 cluster/subspace로 판단하며 arbitrary rotation의 개별 vector를 외부 비교하지 않는다.

## Static partition과 accepted 상태

Linear은 Kff uf=Ff−Kfc uc다. Nonlinear free residual은 fint−F이고 Ktff correction을 풀되
prescribed 값을 factor별 absolute target으로 복구한다. 외력과 prescribed target을 같은 factor로
증가시키며 이전 accepted state로부터 trial을 시작한다. Empty free solve는 유효한 recovery
경로다. Constrained reaction/full residual을 구분하고 모든 required recovery가 끝나야 commit한다.

## HHT constant-mass 식

Public alpha=−.05, beta=(1−alpha)²/4=.275625, gamma=.5−alpha=.55다.
Algorithmic A와 physical endpoint a를 구분한다.

\[
u_p=u_n+hv_n+h^2(.5-\beta)A_n,\quad
v_p=v_n+h(1-\gamma)A_n,
\]
\[
A_{n+1}=(u_{n+1}-u_p)/(\beta h^2),\quad
v_{n+1}=v_p+\gamma hA_{n+1},\quad
a_{n+1}=(A_{n+1}+\alpha a_n)/(1+\alpha).
\]
\[
R=MA_{n+1}+(1+\alpha)(f_{n+1}-F_{n+1})-\alpha(f_n-F_n),
\quad K_{eff}=M/(\beta h^2)+(1+\alpha)K_t.
\]

f/U/Kt는 위 평균 trace energy에서 함께 평가한다. Linear HHT는 그 zero-state K0,
비선형 HHT는 trial의 cross-GP 전체 Hessian이다. M은 상수 reference lumped M_L이며 계수·고정 h·
100회·strict physical/correction/update 판정은 변경하지 않는다.

이전 free endpoint가 \(Ma_n+f_n-F_n=0\)이면 R=(1+alpha)(Ma_next+f_next−F_next)다.
Endpoint 표현을 쓰면 tangent도1/(1+alpha)로 함께 환산해야 한다. Incomplete previous endpoint를
동등하다고 가정하거나 weighted force에서 reaction을 계산하지 않는다. Physical reaction은
Ma+f−F의 constrained projection이다. Initial u/v0에서 Mff a0=Ff(0+)−ff(0), A0=a0다.
곡선은 endpoints에서 실제 평가하며 prescribed dynamic q/v/a/A는 exact0다.

Physical momenta는 current nodal x/v와 같은 reference lumped measure를 사용해
\(p=\sum_a m_av_a\), \(L=\sum_a x_a\times m_av_a\)다.
Rates는 \(\dot p=\sum_a m_aa_a\), \(\dot L=\sum_a x_a\times m_aa_a\)이며
각 node에서 v×v 항은0이다. 독립 oracle도 같은 `rho*N_a*dV0` 적분을 사용한다.
External/support moment는 같은 current 위치로 구한다. HHT 및 force trapezoid work는 frame별
exact energy conservation을 의미하지 않는다. Wext는 endpoint force·du trapezoid 누적,
error=T+U−T0−U0−Wext다. Alpha0 linear unforced benchmark에서는 독립 보존을 확인한다.

Current moment의 검사에는 선형화의 적용 범위를 구분한다. Full physical residual을
\(r=Ma+f_{int}-F\), 공개 constrained reaction을 \(Q=\Pi_c r\), free residual을
\(e=\Pi_f r\)라 하면 \(B_x=\sum x_a\times(F_a+Q_a)-\dot L
=\sum x_a\times f_{int,a}-\sum x_a\times e_a\)다. 정적에서는 관성항을0으로 둔다.
TL은 유한회전 불변성으로 \(\sum x\times f_{int}=0\)지만 선형 K0는
\(\sum X\times K_0u=0\)이며 \(\sum x\times K_0u=\sum u\times K_0u\)가 일반적으로
O(u²)이다. 선형 current 수지 차이를 exact0으로 보정하거나 roundoff gate로 거부하지 않는다.
Current moment와 momentum은 그대로 계산하고 이 차이와 reference-arm identity를 내부
evidence/로그로 설명한다. Newton 허용오차의 free residual 기여도 함께 남긴다. 시간차분
\(\Delta L/h\)를 instantaneous \(\dot L\)와 exact equality로 판정하지 않는다.

## 정식화 검증 경계

Analytic brick/affine energy, independent tensor polynomial integrals,6 rigid modes/rank, scaling/
covariance, node0 포함 gradient/tangent signed directional step sweep와 small-load limit를 결합한다.
Compatible 비균일 trace field에서 reference-weighted mean과 center/unweighted mean의 차이,
cross-GP outer product, 위 weighted-variance energy 차이와 linear–NL zero-state 일치를 검사한다.
GP-local bulk tangent를 사용하는 대조는 방향 미분 검사에서 실패해야 한다.
Warped geometry는 independent higher-order quadrature로 평균과 deviatoric 적분의 discretization
error를 보고하되 production8점 규칙을 바꾸지 않는다. Extra-null·conditioning·굽힘/체적 잠김과
spatial refinement를 별도로 검사한다. Generalized eigen은
작은 independent dense oracle와 original sparse-pencil metrics로 검증한다. HHT는 analytic
oscillator/비선형 constant-mass 제조해에서 u/v/physical a·h/h2/h4 수렴을 검증한다.
외부 comparison의5%/1%로 이 내부 검사를 대신하지 않는다.