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

# Hex8 정식화와 검증 불변량

적용: [PRD](PRD.md)의 FEM-REQ-HEX8-002–007,012–018,021.
물리 선택은 [ADR](ADR.md)에 있다. 이 문서는 승인된 식의 구현용 유도이며
외부 결과에 맞추어 physics를 바꾸는 권한이 아니다. 수치 controls는
[별도 원본](numeric-controls.md), 출력 tensor/identity는 [I/O](io.md)를 따른다.

## 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 재배열로 숨기지 않는다.

## 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를 따른다.

## Linear elasticity와 K0

\(\mu=E/[2(1+\nu)],\lambda=E\nu/[(1+\nu)(1-2\nu)]\).
\(E>0,-1<\nu<.5\)에서 bulk/shear modulus가 양수다. Engineering strain vector는
`[ε11,ε22,ε33,2ε12,2ε13,2ε23]`; stress vector는 실제 shear stress다.
Normal C block은 diagonal λ+2μ/off-diagonal λ, shear diagonal μ다.

\[
K_0=\sum_g B_g^TCB_gw_g,\qquad f=K_0u,\qquad U={1\over2}u^TK_0u.
\]

동일 point의 block 형태는 \(\lambda G_{ai}G_{bk}+\mu\delta_{ik}G_a\cdot G_b+
\mu G_{ak}G_{bi}\)다. B/Voigt 코드와 독립 tensor/polynomial oracle를 구분한다.
Unconstrained K0는 symmetry/PSD 및6 rigid modes를 가지며 BC 적용 후 mechanism/singularity를
solver가 명시적으로 처리한다. Full integration 잠김을 안정화로 바꾸지 않는다.

## Total Lagrangian energy·내력·접선

\[
F=I+\sum_a u_a\otimes G_a,\quad E_G=(F^TF-I)/2,
\quad W={\lambda\over2}(\operatorname{tr}E_G)^2+\mu E_G:E_G.
\]

\(S=\partial W/\partial E_G=\lambda\operatorname{tr}(E_G)I+2\mu E_G\),
\(P=FS\). \(\delta E_G=\operatorname{sym}(F^T\delta F)\)이므로
\(\delta W=P:\delta F\), \(\delta F=\sum\delta u_a\otimes G_a\)다.

\[
U=\sum W_gw_g,\quad f_{ai}=\sum P_{iJ}G_{aJ}w_g.
\]

\(\delta P=\delta F S+F\delta S\)를 분리하면

\[
\mathbb C_{IJLN}=\lambda\delta_{IJ}\delta_{LN}
+\mu(\delta_{IL}\delta_{JN}+\delta_{IN}\delta_{JL}),
\]
\[
(K_{geo})_{ai,bk}=\sum\delta_{ik}G_{aJ}S_{LJ}G_{bL}w_g,
\quad (K_{mat})_{ai,bk}=\sum G_{aJ}F_{iI}\mathbb C_{IJLN}F_{kN}G_{bL}w_g.
\]

Kt=Kmat+Kgeo는 energy Hessian이며 load-independent CLOAD tangent는0이다.
S는 symmetric이고 total Kt도 symmetric이나 stressed trial에서 SPD라고 가정하지 않는다.
Zero state에서 F=I/S=0/geo=0/total=K0다. 유한 rigid R에서 F=R, Green/W/f는0이며
strained 상태를 global R로 회전하면 energy invariant·force/tangent covariance를 만족해야 한다.
Affine deformation의 energy/stress 및 multi-element constant-stress patch가 독립 oracle다.

공식 비선형 stress는 \(\sigma=FSF^T/\det F\), Green strain은 기준 global tensor다.
Linear 공식 strain은 infinitesimal tensor다. 저장 strain의 shear는 실제 tensor 성분이고
B 계산용 engineering shear의 절반이다.8점 값을 그대로 저장한다.

## Consistent mass와 generalized eigen

\[
M_{ai,bj}=\delta_{ij}\sum_g\rho_0N_aN_bw_g,
\quad T={1\over2}v^TMv.
\]

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다. Production에서는 lump하지 않는다.

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와 SPD evidence를 포함해야 하며 이 수학적 식만으로 설치 backend의
실행을 검증했다고 주장하지 않는다.

\(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.
\]

이전 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는 interpolated current x/v를 사용해 \(p=\int\rho v\,dV_0\),
\(L=\int x\times\rho v\,dV_0\)다.8점 consistent quadrature와 nodal mass 표현이 일치해야 한다.
Angular-momentum rate는 \(\int x\times\rho a\,dV_0\)이며 v×v 항은0이다.
Consistent symmetric M에서 nodal 표현의 \(\sum v_a\times(Mv)_a=0\)도 성립한다.
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, scaling/
covariance, gradient/tangent directional step sweep와 small-load limit를 결합한다. Warped geometry는
independent higher-order quadrature로 discretization error를 보고한다. Generalized eigen은
작은 independent dense oracle와 original sparse-pencil metrics로 검증한다. HHT는 analytic
oscillator/비선형 constant-mass 제조해에서 u/v/physical a·h/h2/h4 수렴을 검증한다.
외부 comparison의5%/1%로 이 내부 검사를 대신하지 않는다.