Files
FESADev/docs/hht-implicit-dynamics/numerical-research.md
T

153 lines
15 KiB
Markdown
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
# HHT 수치 제어·검증 설계 추가 근거
- Feature: `hht-implicit-dynamics`
- Owner: `research-agent`
- Date: `2026-09-22`
- Status: `Evidence and verification candidates; numerical thresholds not selected`
- 기준: [프로젝트 지침](../../AGENTS.md), [PRD](../PRD.md), [ADR](../ADR.md#adr-028), [기존 연구](research.md).
[관성 정식화](formulation-inertia.md)와 [독립 probe](numerical-review-probes.md)를 읽고
남은 질문만 조사했다. 기존 SDOF/moving-null 계산은 재실행하지 않았다.
**사실**은 열람한 일차 자료, **도출**은 명시한 수학적 계산, **후보**는 검증 제안이다.
수치 threshold·질량 채택·물리 범위·입출력 계약은 이 문서에서 결정하지 않는다.
현재 질량 선택은 사용자 결정으로 consistent mass가 확정되었다. Shell 속도 보간의
운동에너지를 그대로 적분하며, lumped는 이번 HHT 구현의 지원 요구에서 제외한다.
## 1. Physical acceleration, startup와 DAE 적용 경계
**사실 [N1].** Brüls–Arnold (2008) author postprint pp.3–4, Eq. (5)–(8)은
물리적 `q_ddot`와 보조 acceleration-like variable을 구분하고 보조 초기값을
물리적 초기 가속도와 같게 둔다. 이 문서의 `A`가 논문의 보조 `a`에 대응한다.
`alpha_m=0`, `alpha_f=-alpha`를 대입하면 현재 후보의
`A_(n+1)=(1+alpha)a_(n+1)-alpha a_n`이다. 같은 논문은 variable M일 때 이를
force-level 가중 평형과 동일시할 수 없다고 명시한다. 이 질문에는 기존 출처가 충분하다.
**사실 [N1]: 정리의 가정은 정확히 구분한다.** Printed p.6, “Equivalent model”은
`[[M,-L],[-h_qddot,I-h_y]]` block의 비특이성과 index-3 조건
`Phi_q g*_lambda`의 비특이성을 사용한다. 여기의 `h`는 논문 Eq. (4)의 algebraic
output function이며 시간간격이 아니다. Controller/algebraic coupling을 제거한
특수형에서는 첫 조건이 M의 비특이성으로 줄어든다. 논문 전체가 무조건 “M만
비특이이면 충분”하다고 읽어서는 안 된다. “Initial value condition”의 초기 오차 조건과
“Constrained case”의 정칙성 하에서 얻은 수렴 결과를 singular drilling에 자동 적용할 수 없다.
**도출/수치 제어에의 함의.** 현재 shell의 `psi=N(q)^T g=0`는 논문의 holonomic
`Phi=0`와 동일한 모델임이 아직 증명되지 않았다. 따라서 다음을 분리해서 관찰한다.
| 관찰 대상 | 필요한 검증 후보 | 남은 한계 |
| --- | --- | --- |
| 첫 시점 | `a0`는 우측 물리 평형·초기 대수조건으로 결정한 뒤 `A0=a0`; 첫 step과 이후 오차를 분리 | 작은 Newton residual로 불일관 초기값이 정당화되지 않음 |
| 물리 a와 A | 동일 시각 정확해에 대한 q/v/physical a 차수와 `A-a`를 별도 측정 | A를 물리 출력·각운동량 미분에 대입하면 안 됨 |
| Current null | 매 configuration에서 `MN`, `N^T I`, `psi`, 양의 부분공간 mass 및 `psi_q N` 정칙성 | 이전 시점 null, 고정 UR3, 전체 M의 역행렬로 대체 불가 |
| Hidden constraints | `psi_q v+psi_t`와 `psi_q a+psi_qq[v,v]+2psi_qt v+psi_tt`의 최대·종단 결함과 dt 수렴 | 위치 평형만 만족해도 두 결함이 0인 것은 아님 |
| 반례 coverage | 현재 움직이는 영공간 예제를 full redundant coordinates로 풀어 기존 reduced-angle probe와 대조 | 기존 probe는 위치조건을 정확히 강제했으므로 full-coordinate Newton 구현 검증이 아님 |
| 고정 간격 | 기존 후보에서 dt를 줄여 같은 물리 종단시각 비교 | adaptive dt·startup 재설정은 별도 scope이며 이번 근거로 추가하지 않음 |
**도출.** Chart 내부의 smooth q에서 projector와 null basis는 같은 물리 영공간을 나타내야
한다. Null vector의 부호·정규화 변경 자체와 rank 손실을 구분해야 한다.
`psi_q`에는 `DN`이 들어가므로 basis를 고정한 finite difference는 별도 식을 검증하게 된다.
투영 오차·hidden defect가 큰 결과를 발견하면 정확한 projection/condensation의 채택 여부는
formulation owner에게 돌린다. 여기서는 추가 projection을 자동 제안·채택하지 않는다.
## 2. Double precision, conditioning 및 미분 검증
**사실 [N2–N4].** IEEE double의 정상수 significand는 53 bits다. LAPACK DGECON은
선택한 1/∞ norm의 reciprocal condition number **추정값**을 반환한다. DGESVX는
equilibration 여부를 구별하고, factor pivot이 정확히 0인 경우와 reciprocal condition이
machine precision보다 작은 경우를 서로 다른 INFO로 보고한다. 후자는 해·오차 추정을
계속 계산할 수 있는 warning이다. 이 공개 library 동작은 FESA failure threshold를
자동 승인하지 않으며 vendor routine 채택 요구도 아니다.
**도출.** Binary64에서 1의 다음 표현수 간격 `eps=2^-52`, round-to-nearest 정상 범위의
unit roundoff `u=2^-53`를 구분한다. 문헌의 “machine epsilon” 정의와 library 반환 의미를
확인하지 않고 두 값을 혼용하면 안 된다. 작은 perturbation의 forward error는
대략 `condition × backward error`로 증폭되므로 residual 하나로 해의 정확도를 보장하지 못한다.
| 질문 | 근거에 따른 분리 및 검사 후보 |
| --- | --- |
| Physical mass-null인가? | `n=[0;J_l^-1 d]`, admissible free-space intersection 같은 분석적 운동학을 원본으로 삼는다. `MN`의 수치 결함과 양의 complement의 conditioning을 각각 검사한다. |
| 작은 양의 eigenvalue인가? | 단위·두께 때문에 작은 inertia를 rcond 하나로 drilling null로 바꾸지 않는다. 예상 추가 null/geometry 이상과 positive-space accuracy loss를 구분한다. |
| 어느 rcond인가? | Full singular M의 rcond는 positive inertia 품질을 설명하지 못한다. Positive-space mass, 초기화 algebraic block, effective Newton matrix를 각각 이름 붙이고 norm·scaling을 기록한다. |
| 단위를 바꾸면? | 병진·회전 단위가 섞인 원행렬 condition은 단위 선택에 의존한다. `q=D qhat`, `Rhat=D^T R/Escale`이면 `Khat=D^T K D/Escale`라는 virtual-work 일관 scaling을 대조 후보로 삼는다. 기준 scale의 수치 선택은 NR 소유다. |
| Library equilibration이면 충분한가? | Row/column equilibration은 solver conditioning 도구다. 물리적 nondimensionalization이나 예상 null proof와 동일하지 않다. 원·scaled 행렬 중 무엇을 추정했는지 보존한다. |
**사실 [N5].** GSL의 central differentiation은 truncation과 roundoff를 함께 추정하여
step을 고르며, 5-point와 3-point 차이를 오차 추정에 쓴다. 따라서 가장 작은 step을
항상 정답으로 간주하지 않는 독립 derivative 검사 설계가 공식 수치 library의 접근과 부합한다.
**도출/후보.** Dimensionless direction p에 대해 scaled 잔차의 중심차분과 `Khat p`를
비교하면 일반적인 smooth·비상쇄 구간에서 오차형은 `C_t delta²+C_r u/delta`다.
최적 step의 `u^(1/3)` 크기는 이 가정에서의 점근적 안내일 뿐 공통 고정값이 아니다.
Zero derivative, 큰 residual cancellation, branch 접근에서는 계수와 error floor가 달라진다.
여러 delta를 sweep하여 감소구간·최저 plateau·roundoff 증가구간을 기록하고, 모든 sample이
같은 허용 chart와 geometry 안에 있는지 확인한다. Unit-aware scale을 쓰고 임의
`max(1, norm)`로 작은 문제를 덮지 않는다.
Inertia partial derivative 검사에서는 q/v/physical a 중 하나만 바꾸고 나머지는 고정한다.
Effective tangent 검사에서는 같은 trial q에 연동되는 v, A, physical a를 모두 갱신한다.
두 검사는 서로 다른 미분이므로 섞지 않는다. 작은 dt의 acceleration 재구성에는
위치 차분 오차가 `1/dt²`로 증폭되므로 시간 차수 측정에서 너무 미세한 dt의 plateau를
수학적 차수 저하로 오인하지 않는다. 기준 oracle 정밀도와 nonlinear solve 오차도 분리한다.
## 3. 3차 회전 미분과 운동량·에너지 case 후보
**사실 [N6].** Higham–Relton의 Theorem 3.5/Algorithm 3.6은 고차 Fréchet derivative를
재귀 block matrix function의 우상단 block으로 구한다. Matrix exponential은 smoothness
조건을 만족한다. `P=[phi]_cross` 및 세 방향 행렬에 적용하면 24×24 exponential 하나의
3×3 우상단 block이 3차 혼합 미분을 준다.
**후보.** 이 방식에 `d0`를 곱한 결과를 Rodrigues/회전 Jacobian의 해석 미분과 비교하면
동일 scalar coefficient 도함수를 재사용하는 자기검증을 피할 수 있다. “Exact derivative”는
수학적 항등식을 뜻하며 부동소수점 exp 계산까지 오차가 없다는 뜻은 아니다.
Production 의존성 추가나 24×24 exp의 runtime 사용은 요구하지 않는다.
**자체 도출: 영회전의 유한 다항 oracle.** `H_p=[p]_cross` 등으로 쓰면
\[
D^3d(0)[p,s,r]=\frac16\sum_{\sigma\in S_3}
H_{\sigma(p)}H_{\sigma(s)}H_{\sigma(r)}d_0.
\]
이는 exponential Taylor series의 cubic 항에서 바로 얻는다. 방향 permutation 대칭,
각 방향의 선형성, 동일 방향 3회 `H_p³ d0`를 함께 검사할 수 있다. 비영 phi의
block-exponential oracle와 `D²d`의 방향 중심차분을 추가해 영점 특례만 맞는 오류를 검출한다.
Chart 경계 안쪽, 작은 phi, 서로 비평행한 방향을 포함하되 경계를 넘어 차분하지 않는다.
| ID | 실행 가능한 검증 설계 후보 | 판정에 필요한 한계 |
| --- | --- | --- |
| NR-E01 | 임의 유효 shell 상태에서 spatial rigid velocity를 부여하여 kinetic energy·총 P·총 L을 독립 volume integration과 대조 | Rigid velocity는 순간 운동학 검사다. 회전하는 탄성 shell이 자유운동에서도 강체를 유지한다는 가정은 하지 않음. Consistent/lumped 각각 자신의 measure를 사용 |
| NR-E02 | `I(q,v,a)`를 q/v/a별 미분하고 전체 effective tangent도 별도 directional check | 3차 director 항과 `D(H^T)` 누락을 분리 검출 |
| NR-E03 | 작은 smooth 회전 구간의 자유 shell에 일관된 initial velocity를 부여하거나 짧은 smooth loading 후 unload, P/L/총에너지·work와 dt 수렴 관찰 | 초기조건/시간하중이 CLI 범위 밖이면 kernel/procedure verification 후보로만 유지. 승인 input subset 확대 금지 |
| NR-E04 | 같은 spatial problem에 alpha=0 및 승인된 음수 alpha를 적용하고 저주파 해상 case와 높은 modal content case를 분리 | 일반 비선형 HHT의 총에너지 매 step 감소 또는 정확한 각운동량 보존을 무조건 요구하지 않음 |
| NR-E05 | 문헌 free-flying single-layer plate를 후속 benchmark 후보로 검토 | 아래 N7은 large rotation/energy 사례의 출처다. 현재 fixed chart를 넘어서는 tumble·장시간 운동은 그대로 채택할 수 없음 |
**사실 [N7].** Vu-Quoc–Tan (2003), §5.1.2/Figs.9–12는 free-flying single-layer plate의
큰 전체 운동과 P/L/에너지 및 dissipative/nondissipative integrator 결과를 비교한다.
해당 요소는 rotation DOF가 없는 solid-shell이므로 그 mass·응력 모델과 여기의 director shell은
다르다. Author-uploaded 본문의 관련 절을 확인했으나 그림의 하중 위치·수치 history를
독립 추출하지 않았다. 따라서 현 단계에서는 case family 근거이며 완성된 reference artifact가 아니다.
**추론.** 현재 shell 운동학에서는 `Pdot/Ldot`에 physical a를 사용하고 numerical drill
potential을 포함한 에너지·내력 일을 일관되게 대조해야 한다. N7의 보존형 integrator가
보이는 exact conservation을 HHT의 요구조건으로 옮길 수 없다. [기존 연구 S8](research.md#sources)의
비선형 Newmark 보존 한계 근거는 그대로 충분하다. 추가 shell 연구는 기하·하중·시간 구간이
승인된 chart와 입력 범위 안에 있는지 결정한 뒤에만 구체화할 필요가 있다.
## 출처·열람 한계와 인계
`accessed: 2026-09-22`. 문서 정적 검사만 수행하며 solver/build/test/reference 생성은 수행하지 않았다.
| ID | Stable citation / locator | 열람 범위 |
| --- | --- | --- |
| N1 | Brüls, O.; Arnold, M. (2008), *The generalized-alpha scheme as a linear multistep integrator: Toward a general mechatronic simulator*, J. Computational and Nonlinear Dynamics 3, 041007, [DOI:10.1115/1.2960475](https://doi.org/10.1115/1.2960475), [author PDF](https://orbi.uliege.be/bitstream/2268/27088/1/pa_BrulsArnold2008_TheGeneralizedAlphaSchemeAsALinearMultistepIntegrator.pdf) | Printed pp.3–6 Eq.5–8/Equivalent model, Initial value condition/Constrained case. Shell singular-mass 정리 아님 |
| N2 | Oracle, *Numerical Computation Guide*, [§2.2.3 Double Format](https://docs.oracle.com/cd/E37069_01/html/E39019/z4000ac019269.html) | 53-bit significand 설명. IEEE 표준 전체 적합성 시험을 수행한 것은 아님 |
| N3 | LAPACK, [DGECON](https://netlib.org/lapack/explore-html/d4/daf/group__gecon_ga4f9b830e19e12c7f082ddb497a57af18.html), Purpose/NORM/ANORM/RCOND | Normwise reciprocal condition estimate |
| N4 | LAPACK, [DGESVX](https://netlib.org/lapack/explore-html/d5/dbe/group__gesvx_ga82173a93234afc15d70b64233b3e5bc8.html), Description steps 1–6/INFO | Equilibration, singular pivot와 working-precision warning 구별 |
| N5 | GNU Scientific Library 2.8, [Numerical Differentiation](https://www.gnu.org/software/gsl/doc/html/diff.html), `gsl_deriv_central` | Truncation/roundoff 기반 step/error estimate. FESA tolerance 원본 아님 |
| N6 | Higham, N.J.; Relton, S.D. (2014), *Higher Order Fréchet Derivatives of Matrix Functions and the Level-2 Condition Number*, SIAM J. Matrix Anal. Appl. 35, 1019–1037, [DOI:10.1137/130945259](https://doi.org/10.1137/130945259), [author PDF](https://eprints.maths.manchester.ac.uk/2066/1/paper.pdf), Theorem 3.5/Algorithm 3.6 | 2013 author version printed pp.5–6; oracle identity와 계산법 |
| N7 | Vu-Quoc, L.; Tan, X.G. (2003), *Optimal solid shells for non-linear analyses of multilayer composites. II. Dynamics*, CMAME 192, 1017–1059, [DOI:10.1016/S0045-7825(02)00336-5](https://doi.org/10.1016/S0045-7825(02)00336-5), [author-uploaded text](https://www.researchgate.net/publication/222546721_Optimal_solid_shells_for_non-linear_analyses_of_multilayer_composites_II_Dynamics), §5.1.2/Figs.9–12 | Relevant text 확인; 하중 figure·곡선 digitization 미수행, exact benchmark reproduction 미정 |
Numerical Review에 인계할 결정은 scaling·conditioning/hidden-defect 판정·FD sweep과
시간 수렴 구간·case 최소 조합이다. Formulation에 남는 사항은 singular-mass regularity의
실제 shell 적용과 full-coordinate hidden constraint 검증이다. 사용자 결정에 따라 이번
HHT 구현은 consistent mass를 사용한다. 위 lumped 비교 근거는 연구 이력이며 현재 지원 요구가 아니다.