docs: consolidate solver knowledge and preserve formulation originals
This commit is contained in:
KOKO\Mimi
@@ -1,70 +1,129 @@
|
||||
# HHT implicit dynamics
|
||||
|
||||
## 목적과 현재 상태
|
||||
## 범위와 현재 상태
|
||||
|
||||
2026-10-11 [HHT-ADR-012](ADR.md)에 따라 nodal lumped kinetic mass로 변경한다.
|
||||
변경본의 관성·convective 항·접선·state/HDF5 관련 Debug 검사38개와 Python 비교기137개가
|
||||
통과했고, 같은 production 소스의 Release 원본 U 101-frame 비교 PASS·UR 경고200건을
|
||||
확인했다. 독립 수치·물리 검토는 bounded PASS다. [채택 보고서](../hex8-structural-analysis/lumped-mass-adoption.md)에
|
||||
검증 구성과 한계를 기록한다. 새 전체 Debug CTest789개는 자연 종료했다(native8:
|
||||
784 PASS·4 FAIL·1 SKIP). Stale Hex8 test 두 개는 수정 후 공식 focused 검사에서 PASS다.
|
||||
이 실행 당시 C08a 시간 정확도 FAIL·C08b 새 Newton 회전 예외 및 C04a SKIP이 남았다.
|
||||
2026-10-12 사용자 요청으로 [두 시작점 정책](newton-start-policy.md)을 구현했다. 관련37개
|
||||
Debug 회귀와 독립 수치 검토가 완료됐다. 전체802 CTest는 제외 없이 자연 종료했다
|
||||
(11384.61초·native8:799 PASS·C08 정확도 FAIL2·기존 C04a SKIP1). 실제 C08b의5개 grid와
|
||||
새65-frame runtime 회귀에서 회전 예외가 해소됐으며 fresh 원본 U101 families PASS·
|
||||
UR 경고200건과 독립 kinetic/balance 대조를 확인했다. Source·원본·binary identity는
|
||||
검증 전후 일치한다. Runtime 해결을 시간 정확도 해결로 해석하지 않는다.
|
||||
다음 기존 consistent 구현의 이력과 원본 변위 통과를 현재 전체 검증으로 승계하지 않는다. 실제 shell의
|
||||
C08a/C08b 시간수렴 실패는 미해결이다. 구현됨과 기능 검증 완료를 구분한다.
|
||||
MITC4 shell의 유한회전 기하비선형 동해석이다. S4/S4R은 FESA-MITC4로 해석하며
|
||||
외부 solver 내부 알고리즘과의 동등성을 뜻하지 않는다. 현재 채택 경로는
|
||||
**nodal-lumped-kinetic 질량, reduced-coordinate HHT, 초기 null의 고정 직교 보완기저,
|
||||
shell schema v2**다. 과거 consistent 질량·full-coordinate v1 근거를 변경본 통과로 승계하지 않는다.
|
||||
|
||||
## 계약과 사용법
|
||||
2026-10-13 기록에서 ADR-018의 실제 C09와 공식 Release C08a 7-grid는 통과했다.
|
||||
최신 전체 Debug는 856개 중 44개 PASS 후 사용자 중단됐고 Release C08b도 중단됐다.
|
||||
당시 원본 Release 해석은 증분99·t=9.9까지 accepted 로그만 남아 HDF5 비교가 불가능했다.
|
||||
사용자 중단 exit -1은 수치 FAIL이나 PASS가 아니다. 이후 2026-10-14 성능 baseline은
|
||||
원본을6372.325621초/native0으로 완료하고 schema v2·729nodes/101frames를 확인했다.
|
||||
9-node U blocking 비교는 PASS(최대 relative RMS .0022824245), UR은 경고200건이었다.
|
||||
이는 baseline source의 결과다. 최신 최적화 source의 새 원본 비교·C08b·전체 회귀·독립
|
||||
전체 물리 판정은 남아 있으며 전체 green이 아니다. [성능 요약](../calculation-performance/README.md)을 함께 본다.
|
||||
|
||||
2026-10-12 후속 [HHT-ADR-014](ADR.md)로 P 판정·C08b 최종 finer grid와 reduced shell
|
||||
HHT/schema v2 구현을 승인받았다. 실제 채택식은
|
||||
[축소 좌표 정식화](formulation-reduced-coordinate.md),
|
||||
[carried-A arithmetic/noise](formulation-carried-a-noise.md),
|
||||
[schema v2 상세](io-reduced-coordinate.md)를 따른다. 현재 구현과 검증은 진행 중이며
|
||||
기존 v1 결과·실패 이력을 변경본 통과로 승계하지 않는다.
|
||||
## 채택 방법과 산술 의미
|
||||
|
||||
- [PRD](PRD.md)·[ADR](ADR.md)·[요구사항](requirements.md): 범위·미해결·현재 checklist.
|
||||
- [정식화](formulation.md), [유한회전 관성](formulation-inertia.md),
|
||||
[가속도 복구 후보](formulation-acceleration-recovery.md): 채택식과 미채택 후보 구분.
|
||||
- [I/O](io-contract.md), [reference 계약](reference-model.md),
|
||||
[전체 변위 비교표](displacement-comparison-tables.md).
|
||||
- [수치 기준](numerical-controls.md), [검증 사례](verification-cases.md),
|
||||
[수학 검산](numerical-review.md#numerical-review), [연구](research.md#research).
|
||||
- [구현 검토·실패·수정](implementation-numerical-review.md), [물리 평가](physics-evaluation.md).
|
||||
- [Consistent/lumped 질량 비교 실험](mass-comparison.md): 채택 전 private 비교 이력과 후속 승인; production 재검증과 구별.
|
||||
- [Newton 회전 실패 원인·해결 후보](newton-rotation-investigation.md): 변경 전 moving-null 접선 상쇄의 독립 대조와 private 후보 실험.
|
||||
- [두 시작점 정책의 구현·검증](newton-start-policy.md): HHT-ADR-013 실제 정책·시행착오·현재 evidence·적용 한계. C08 시간 정확도와 구분.
|
||||
- [시간 정확도·P 감소 실패의 해결 후보 조사](research.md#temporal-remediation-research): 축소 좌표·속도 수준 DAE 적분·시간 해상도·합산의 조사와 채택 경계. Solver 변경 또는 실패 폐쇄를 뜻하지 않는다.
|
||||
- [시간 정확도·P 후속 검증](temporal-remediation.md): 소형 DAE 독립 검산, 실제 finer-grid와 carried-A private 산술 후보, P strict-decrease 계약 문제. 기존 C08 FAIL과 미채택 알고리즘을 구분한다.
|
||||
- [Reduced C08 구현](implementation-c08-reduced.md): 실제 controller·N7·전체 산술 이력 및 좌표 증분 누적의 경계.
|
||||
- [C08a 추가 refinement 결정](c08a-refinement-proposal.md): HHT-ADR-015로 승인된 두 변경과 승인 전 진단 evidence. 변경본 공식 검증은 별도로 수행한다.
|
||||
- [N3 좌표 표현 결정](c08-n3-coordinate-representation-proposal.md): HHT-ADR-016으로 승인된 exact pair 증분 검사와 RN-even 저장 검사. 승인과 변경본 검증 완료를 구분한다.
|
||||
- [병진 유지 Newton 명세](c09-channel-newton-proposal.md), [독립 검토와 N7 대응](numerical-review-channel-newton.md): HHT-ADR-017로 승인된 제한 보정·부분 행렬 N2-L 검사·full-model defect 보존. [C09 산술 진단](c09-newton-rounding-diagnostic.md)의 endpoint 존재증거와 자연 trajectory 검증을 구분한다.
|
||||
- [Channel Newton 구현과 실행 기록](implementation-channel-newton.md): private helper/자연 proper-Q 진단과 configured 검증의 범위를 구분한다.
|
||||
- [유한회전 공변성 실패 진단](channel-newton-covariance-diagnostic.md), [초기 null 직교 보완공간 제안](fixed-null-basis-proposal.md): 실제 selector recurrence 반례와 HHT-ADR-018로 승인된 basis/integrator 변경의 경계.
|
||||
- [고정 null 기저 상세 정식화](formulation-fixed-null-basis.md), [독립 수치 검토](numerical-review-fixed-null-basis.md): sparse projection/lift, 혼합 부분구속 chart, augmented scaled operator, 초기 상태와 CA/N7의 구현 계약. 정적 수식 검토와 변경본 실행 검증은 별도다.
|
||||
- [고정 null 기저 구현과 실행 기록](implementation-fixed-null-basis.md): 실제 scaled solve, basis 계수 오차 반경, global 좌표 반복 감지와 현재 검증의 범위.
|
||||
- [고정 null 기저 물리 평가](physics-evaluation-fixed-null-basis.md): 새 C09의 실제 검증 범위, 강제 강체 병진의 독립 한계 계산과 reference 물리 검산의 진행 상태.
|
||||
- [사용자 중단 시점 비교](interrupted-validation-comparison.md): Debug44개 완료와 Release 요약 일치, 원본99증분의 수렴 로그, 현재 HDF5 부재로 인한 Abaqus 비교 한계 및 후속 실행 순서.
|
||||
[기본 정식화](formulation.md), [관성](formulation-inertia.md),
|
||||
[축소 좌표](formulation-reduced-coordinate.md), [고정 null 기저](formulation-fixed-null-basis.md)가
|
||||
수식의 상세 원본이다. 초기 null과 실제 free 공간의 교집합에서 고정 sparse 열 Lp/Lz를
|
||||
만들고 coupled graph root를 풀어 full `q=Q(p,t)`를 재구성한다. v와 physical a는 Q의
|
||||
chain rule로 복원한다. 부분 회전구속을 full clamp로 바꾸거나 모든 node에 동일 reduced
|
||||
차원을 강제하지 않는다. Fixed patch 밖에서 검증된 history transport가 없으면 실패한다.
|
||||
|
||||
환경·명령은 [operations](../operations.md)를 따른다. I/O의 지원 `*DYNAMIC, DIRECT` subset과 고정 HHT 계수,
|
||||
출력의 물리적 가속도와 보조 acceleration 구분을 확인한다.
|
||||
HHT 계수는 `alpha=-.05,beta=.275625,gamma=.55`다. Private p/vp/ap/Ar가 적분 이력이며
|
||||
`Ar_n=(1+alpha)ap_n-alpha*ap_(n-1)`다. 공개 full `A_n=(1+alpha)a_n-alpha*a_(n-1)`는
|
||||
**full-chart 진단용 filter combination**이고 full q/v의 predictor driver가 아니다.
|
||||
Frame0은 `q0=v0=0,A0=a0`이며 full 초기평형·hidden consistency를 만족해야 한다.
|
||||
공개 A 또는 Qp*Ar로 private history나 restart를 구성할 수 없다.
|
||||
|
||||
## 재사용 Handoff
|
||||
관성력 `I=M(q)a+c(q,v)`, 접선, 에너지·운동량은 같은 nodal kinetic measure에서 계산한다.
|
||||
Lumped는 회전 block의 scalar diagonal 치환이 아니다. Epsilon inertia로 null을 없애지 않는다.
|
||||
Primitive hidden-gradient authority, tiny director increment와 coefficient-first inertia contraction은
|
||||
큰 병진과 작은 회전의 소실을 줄이는 구현이다. Force·tangent·recovery는 같은 state/energy를
|
||||
사용한다. [Carried-A 산술](formulation-carried-a-noise.md)은 noise indicator의 범위를,
|
||||
[가속도 복구 후보](formulation-acceleration-recovery.md)는 채택 경로와 다른 후보를 설명한다.
|
||||
|
||||
Eigen의 기존 nodal block lumped mass는 유지하고 HHT도 동일한 scalar moments의 nodal kinetic measure를 사용한다. HHT의 회전 의존 M·convective force·접선과 운동량은 같은 measure로 계산하며 scalar diagonal 치환은 하지 않는다. Total rotation-vector
|
||||
chart의 generalized moment와 spatial moment도 다르다. Physical a와 algorithmic A를 바꾸어
|
||||
쓰지 않으며 moving-null 제약과 동적·대수 초기평형을 유지한다. 같은 trial에서 full residual,
|
||||
inertia와 tangent를 계산하고 실패 candidate가 accepted history를 변경하지 않게 한다.
|
||||
## 입력·출력·reference 계약
|
||||
|
||||
## 실패·시행착오·다음 작업
|
||||
Single Part/identity Instance, homogeneous centered shell, 상수 isotropic elastic material과
|
||||
양의 DENSITY, NLGEOM=YES인 단일 `*DYNAMIC, DIRECT`를 지원한다. INC 기본값은100이며
|
||||
양의 uint32 증분 상한이다. h<=T이고 double(N)*h==T, N<=INC를 요구한다. h/T/INC는 정수 증분 수와 binary64
|
||||
n*h의 exact schedule로 검증한다. 영 BC·영 초기 q/v·일정 global 병진 CLOAD만 허용한다.
|
||||
사용자 계수, adaptive/explicit, amplitude/follower/분포하중, 회전 CLOAD, 비영/시간변화 BC,
|
||||
비영 초기조건, restart는 제외한다. 허용된 `*RESTART, WRITE, FREQUENCY=0`은 warning no-op이다.
|
||||
CLI·lease·공통 identity/비교 정책은 [solver 공통 요약](../solver-core/README.md)을 따른다.
|
||||
|
||||
C08 두 실패는 열해석 전후 전체617개 회귀에서도 baseline metric과 동일했다.
|
||||
가속도 복구 후보는 아직 채택하지 않았다. C09–C11의 잔여 coverage와 전체 physics 완료도
|
||||
남아 있다. 굽힘 단면력의 두께 계수 `M=t^3 C kappa/12` 수정은14 focused tests 및 원본
|
||||
HDF5 대조의 별도 evidence이며 동해석 전체의 완료를 뜻하지 않는다.
|
||||
과거475 PASS/1 SKIP의 제외 실행과 최신 전체 회귀를 합치지 않는다.
|
||||
Schema identity는 `(hht-implicit-dynamics,hht-implicit-dynamics,2)`다.
|
||||
`/metadata::time_integration_formulation`은
|
||||
`reduced-coordinate-HHT-with-physical-rate-reconstruction`,
|
||||
`mass_formulation=nodal-lumped-kinetic`,
|
||||
`acceleration_convention=physical-endpoint-and-auxiliary-separated`다.
|
||||
모든 result dataset version도 2이며 A의 quantity_semantics는
|
||||
`full-chart-HHT-filter-combination-not-predictor-driver`다. 잘못된 marker·혼합 version·mass·
|
||||
inventory·nonfinite는 dynamic-output-schema 실패다. Historical v1 reader branch는
|
||||
과거 full-coordinate 결과만 읽으며 marker를 보충하거나 v2로 재표시하지 않는다.
|
||||
|
||||
각 `/steps/Step-1/frames/<n>`의 주요 F64 little-endian 결과는 다음과 같다.
|
||||
|
||||
| 상대 경로 | Shape와 의미 |
|
||||
| --- | --- |
|
||||
| nodal/displacement, velocity, acceleration, algorithmic_acceleration | [N,6]; translation/global total rotation-vector chart의 q/v/physical a/public A |
|
||||
| nodal/angular_velocity, angular_acceleration | [N,3]; Jl*phidot 및 Jl*phiddot+DJl[phidot]phidot |
|
||||
| nodal/triad, director | [N,3,3], [N,3] |
|
||||
| nodal/generalized_inertia_force, generalized_internal_force, generalized_external_force, generalized_residual, generalized_reaction | [N,6]; residual=I+f-F, generalized reaction=constrained projection |
|
||||
| nodal/reaction | [N,6]; support-on-structure spatial RF/RM, moment는 Jl^-T 변환 |
|
||||
| element/shell/local_frame, generalized_strain, section_resultant, stress | [E,4,3,3], [E,4,8], [E,4,8], [E,4,3,3]; stress는 bottom/middle/top × S11/S22/S12 |
|
||||
| global/energy, global/equilibrium | [5] kinetic/physical/drill/work/defect, [6] 관성 포함 force/moment balance |
|
||||
|
||||
Row는 model source inventory와 결합한다. Generalized moment와 spatial moment, physical a와
|
||||
A, drilling force와 physical stress를 혼합하지 않는다. 영 BC는 q/v/a/A exact zero다.
|
||||
일정 외력의 work는 `F*^T(qn-q0)`, defect는 `T+Uphysical+Udrill-E0-W`다.
|
||||
Momentum·hidden 진단 전체가 HDF5 필수 dataset인 것은 아니다.
|
||||
|
||||
Reference는 reference/shellNL-dyn-meshfine의 INP와 displacements CSV다. 전체729-node/
|
||||
101-frame identity·finite를 먼저 확인한 뒤 고정9-node subset으로 투영한다. Instance/node/
|
||||
Step-1/increment0..100으로 join하며 CSV decimal n/10과 binary double(n)*double(.1)을 각각
|
||||
검증한다. U는 displacement 앞3성분, UR은 뒤3성분이다. Frame별27 scalars의 U blocking·
|
||||
UR warning-only 총202 families다. S=max|reference|, near-zero band/error limit .01S,
|
||||
일반 상대오차 .05, family RMS/S .01, S=0 exact-zero다. Warning-only여도 schema/identity/
|
||||
nonfinite 오류는 blocking이다. CSV에는 V/a/A/energy가 없어 변위 일치로 이들을 검증할 수 없다.
|
||||
|
||||
## 코드 흐름과 수명
|
||||
|
||||
[HhtDynamicProcedureExecutor::Execute](../../src/fesa/analysis/hht_dynamic_procedure.cpp)는
|
||||
[ExecuteReducedHht](../../src/fesa/analysis/hht_reduced_procedure.cpp)로 연결된다.
|
||||
[FixedChart/FixedBasisColumn](../../src/fesa/analysis/hht_reduced_internal.h)과
|
||||
[BuildFixedChart·SolveEmbeddingRoot·rate/augmented operators](../../src/fesa/analysis/hht_reduced_operators.cpp)는
|
||||
analysis-private다. Sparse support가 projection/lift authority이며 DOF label만으로 단일
|
||||
성분 injection을 가정하지 않는다. Domain은 non-owning AnalysisModel보다 오래 살아야 하고
|
||||
호출자 소유 solver/parallel adapter는 Execute 동안 유효해야 한다.
|
||||
|
||||
Private chart/history와 [DynamicAnalysisState](../../include/fesa/analysis/dynamic_analysis_state.h)의
|
||||
full accepted frame은 같은 acceptance 뒤 commit한다. Full residual·recovery·finite·constraint
|
||||
검증 전에 일부 이력을 commit하지 않는다. [DynamicHdf5ResultsWriter](../../src/fesa/io/hdf5/dynamic_hdf5_results_writer.cpp)는
|
||||
완결 이력을 write/close/reopen/self-check한 뒤 final을 교체한다. 실패 candidate는 기존 final을
|
||||
보존한다. Accepted history는 실행 중 메모리에 있어 강제 종료 로그만으로 복구되지 않는다.
|
||||
|
||||
## 검증 기준·교훈·남은 한계
|
||||
|
||||
N2 full force/moment·q/v/a/A correction은 채널별 1e-10, 실제 scaled solve backward error는
|
||||
1e-12다. N3 update는 상쇄 전 same-unit scale의 1e-12와 exact-zero branch를 유지한다.
|
||||
Pair increment의 exact expansion 검사와 RN-even 저장 검사는 별도다. Hidden 위치 조건·rank/
|
||||
regularity를 유지하되 hidden rate 진단을 임의 Newton gate로 바꾸지 않는다. 일반 dynamic
|
||||
operator에 대칭성을 강제하지 않는다. 일반 nonlinear HHT energy에 매 frame 보존·감소를
|
||||
요구하거나 numerical dissipation을 physical damping으로 기록하지 않는다.
|
||||
|
||||
시간 검증은 q/v/physical a의 마지막 두 refinement 감소·차수 >=1.8 및 finest normalized
|
||||
error <=1e-4를 요구한다. N7 read-only linear/nonlinear defect solves와 arithmetic indicator는
|
||||
100배 separation을 확인하며 state polishing이 아니다. 원본의 성긴 dt/T1≈7.4567에서
|
||||
변위 일치와 시간 정확도를 구분한다.
|
||||
|
||||
과거 전체617의 HHT2 FAIL/1 SKIP, lumped 전체789의784 PASS/4 FAIL/1 SKIP, 후속802의
|
||||
799 PASS/2 FAIL/1 SKIP은 각 당시 source 이력이다. 초기 C08b runtime 예외는 두 시작점
|
||||
정책으로 해소됐고 후속 C09/C08a7 통과와 구분한다. 최신 Debug h=1/64 완료 요약의 이전
|
||||
Release 일치는 전체 raw state bitwise나 전체 시간수렴 증거가 아니다. 원본 accepted residual
|
||||
최대 병진1.889238010e-15·회전1.452384011e-15도 Abaqus 변위 오차가 아니다.
|
||||
C08b finer-grid·최신 최적화 source의 원본 HDF5 비교·전체 회귀·전체 물리 closure가 남는다.
|
||||
|
||||
## 보존한 과거 수식과 연구의 위치
|
||||
|
||||
| 원본 | 현재 해석 |
|
||||
| --- | --- |
|
||||
| [2026-10-04 본식](formulation-historical-docs-cleanup-20261004.md), [관성 원본](formulation-inertia-historical-docs-cleanup-20261004.md) | 당시 consistent 정책의 수식·검토 이력이다. 현재 lumped/reduced/fixed-null 정책이 우선한다. |
|
||||
| [C08b Newton 연구](formulation-c08b-newton-research-20261012.md), [DAE 연구](formulation-dae-research-20261012.md) | 원인·후보 조사이며 모든 Newton/DAE 후보의 채택을 뜻하지 않는다. |
|
||||
| [운동량 연구](formulation-momentum-research-20261012.md), [graph certificate 초안](formulation-graph-certificate-draft-20261012.md) | 불변량·controller/certificate 연구의 원본이다. 현재 acceptance나 엄밀한 전역 오차 보증을 자동 대체하지 않는다. |
|
||||
@@ -0,0 +1,118 @@
|
||||
# C08b Newton 잔차 재정식화 후보
|
||||
|
||||
상태: 연구 후보이며 production 채택·수치 정확도 승인이 아니다. 2026-10-12 formulation owner 작업은 이 문서만 소유한다. 질량·시간적분·내력·초기화·tolerance·원본을 바꾸지 않는다. 독립 numerical review와 실제 prototype 결과는 별도 owner가 소유한다.
|
||||
|
||||
계약: [기존 정식화](../../docs/hht-implicit-dynamics/formulation.md), [관성·null 상세 §4](../../docs/hht-implicit-dynamics/formulation-inertia.md#4-mass-null과-초기-일관성), [ADR-010/012](../../docs/hht-implicit-dynamics/ADR.md). 실제 trial은 [analysis.json](analysis.json)에 있으며 아래 기전 설명은 그 수치를 바탕으로 한 조사 후보이다.
|
||||
|
||||
## 1. 적용 가정과 보존량
|
||||
|
||||
우선 대상은 9 nodes, 54 DOFs가 모두 free인 C08b 평판 fixture이다. 절점 순서는 기존 source identity이며 각 node의 좌표는 `q_i=(u_x,u_y,u_z,phi_x,phi_y,phi_z)`다. 시간별 미지수는 total rotation chart의 `q`; `v(q),a(q),A(q)`는 기존 HHT/Newmark 종속량이다. Director는 `d=exp([phi]x)d0`, `B=D_phi d=-[d]x J_l`이며 `|phi|<pi`와 모든 기존 geometry 조건을 유지한다.
|
||||
|
||||
채택한 nodal kinetic mass는 rank 5/node이며 inertia covector는 `I=M(q)a+c(q,v)`다. Physical/drilling energy, constitutive law, quadrature, force sign 및 recovery는 그대로 유지한다. 아래 식의 `K_int`는 physical와 numerical drilling을 모두 포함한다. C08b의 force는 0이며 follower/시간변화 force가 없다.
|
||||
|
||||
Endpoint 방정식과 정확한 total derivative를 정의한다.
|
||||
|
||||
\[
|
||||
E(q)=I(q,v(q),a(q))+g(q),\qquad g=f_{int}-F,
|
||||
\]
|
||||
\[
|
||||
K_{eff}=D_qE=K_{int}+D_qI+\frac{\gamma}{\beta h}D_vI+
|
||||
\frac{M}{(1+\alpha)\beta h^2}.
|
||||
\]
|
||||
|
||||
`E_u`는 force, `E_phi`는 force·length이다. `phi`는 rad, `v_phi,a_phi`는 각각 1/time, 1/time²이다. 새 equation은 회전 residual의 단위를 바꾸지 않는다.
|
||||
|
||||
## 2. Moving null 미분과 raw Newton 약화의 후보 기전
|
||||
|
||||
각 node에 대해 rotational null과 단위 null을 정의한다.
|
||||
|
||||
\[
|
||||
\nu_i=J_l(\phi_i)^{-1}d_i,\qquad
|
||||
n_i=\nu_i/|\nu_i|,\qquad
|
||||
N_i=[0;n_i]
|
||||
\]
|
||||
|
||||
`N`의 열은 서로 다른 node에만 지지하므로 all-free fixture에서는 orthonormal이다. `B_i nu_i=0`, `M N=0`, `N^T I=0`는 같은 state의 kinetic mapping에서 성립한다. 임의 방향 `p`에 대한 analytic derivative는
|
||||
|
||||
\[
|
||||
D\nu_i[p]=J_l^{-1}\{B_i p_{\phi i}-DJ_l[p_{\phi i}]\nu_i\},
|
||||
\]
|
||||
\[
|
||||
Dn_i[p]=\frac{(I_3-n_in_i^T)D\nu_i[p]}{|\nu_i|}.
|
||||
\]
|
||||
|
||||
Translation 방향은 `Dn=0`이다. 작은 각에서는 기존 analytic series를 사용한다. Null 부호를 매 평가마다 임의 뒤집거나 unnormalized null의 미분을 unit null에 그대로 사용하지 않는다. `psi` 자체를 다시 scaling/normalization하면 그 configuration derivative도 필요하다.
|
||||
|
||||
연속 대수 residual은 `psi=N^T g`이고 `N^T E=psi`이다. 이 identity를 HHT trial mapping을 포함해 미분하면
|
||||
|
||||
\[
|
||||
\boxed{N^T K_{eff}=D\psi-(DN)^TE.}
|
||||
\]
|
||||
|
||||
따라서 큰 off-root `E`에서는 `-(DN)^T E`가 실제 algebraic derivative `Dpsi`를 상쇄할 수 있다. 이것은 정확한 raw Jacobian 자체에서도 가능한 Newton conditioning 기전이며, 위 identity가 검증됐다고 kernel bug가 증명되는 것은 아니다.
|
||||
|
||||
Newton 시작을 `q_n`으로 두면 `q_n-q_p=O(h)`이며 `A=(q-q_p)/(beta h²)`는 일반적으로 `O(1/h)`인 시험 가속도를 만든다. 이것은 accepted acceleration이 발산했다는 뜻이 아니다. 작은 h에서 off-root inertia residual과 moving-null connection 항이 커질 수 있는 한 경로다. 현재 inc56의 큰 null correction은 이 경로와 일치하는 후보 evidence이며 최종 인과 판정은 독립 numerical owner가 수행한다.
|
||||
|
||||
## 3. 같은 endpoint root를 갖는 정확한 잔차 후보
|
||||
|
||||
`P`는 각 node의 `ux,uy,uz,phi_x,phi_y` equation을 선택하는 고정 45×54 matrix이다. 이는 coordinate row 선택이지 positive-mass 공간으로의 직교 projector가 아니다. 다음 54 equations를 푼다.
|
||||
|
||||
\[
|
||||
\boxed{\mathcal F(q)=\begin{bmatrix}PE(q)\\\psi(q)\end{bmatrix},\qquad
|
||||
\mathcal J(q)=\begin{bmatrix}P K_{eff}(q)\\D\psi(q)\end{bmatrix}.}
|
||||
\]
|
||||
\[
|
||||
\boxed{D\psi[p]=(DN[p])^Tg+N^TK_{int}p.}
|
||||
\]
|
||||
|
||||
현재 고정 force에서는 위 식이 정확하다. Configuration-dependent external force를 지원하는 다른 계약에서는 `K_int`를 `Dg=K_int-D_qF`로 바꿔야 하며 현재 후보의 지원 범위를 자동 확대하지 않는다. Prototype의 finite difference derivative는 위 analytic 식의 독립 검증용이며 production derivative 채택이 아니다.
|
||||
|
||||
`T(q)=[P;N(q)^T]`라 두면 `mathcal F=T E`이다. 현재 선택에서는 각 `n_i,z != 0`가 필요하고 충분한 가역 조건이다. 유효한 지역에서
|
||||
|
||||
\[
|
||||
\mathcal F=0\iff E=0.
|
||||
\]
|
||||
|
||||
보존되는 것은 정확한 root 집합이다. Root 근처에는 `mathcal J=T K_eff`; off-root에는 null rows에 `(DN)^T E`가 추가되므로 Newton correction은 달라진다. 단순한 row preconditioning `T K_eff delta=-T E`는 가역 T 아래 raw correction과 같아 이 후보 효과를 만들지 못한다.
|
||||
|
||||
후보는 `mathcal J delta q=-mathcal F`, `q<-q+delta q`의 full Newton이다. 차분 tangent나 frozen null 미분을 생략한 modified Newton으로 해석하지 않는다. 대수 row는 artificial mass·drill torque·물리 damping을 추가하지 않는다. 다만 기존 `E`에 대한 full Newton에서 다른 nonlinear residual에 대한 full Newton으로 iteration 정책이 바뀌므로 ADR-010의 명시 계약과 별도 채택 경계가 필요하다.
|
||||
|
||||
가역성이 충분히 유지되고 같은 지역 root로 수렴한 경우 discrete endpoint/history는 기존 방정식의 해다. Multiple roots와 branch selection이 있는 비선형 문제에서는 root 집합 동일성이 같은 branch 선택까지 보장하지 않는다.
|
||||
|
||||
## 4. Scaling·판정·state update
|
||||
|
||||
정확한 arithmetic에서 nonsingular constant row/column scaling은 correction을 바꾸지 않지만 floating-point conditioning은 달라진다. Translation과 rotation의 단위가 다르므로 raw unscaled condition number의 절댓값을 보편적 물리 기준으로 사용하지 않는다. Existing scale-aware 검사와 numerical review의 dimensionless 분석을 병행한다.
|
||||
|
||||
매 linear solve에서 `W mathcal J C y=-W mathcal F`, `delta q=C y`와 같이 같은 W/C를 적용하는 algebraic scaling은 허용 연구이다. 이를 `W(q) mathcal F(q)`라는 새 nonlinear residual로 정의하면 `DW·mathcal F`가 정확한 derivative에 필요하다. 두 해석을 혼합하지 않는다. Null normalization derivative는 위 §2처럼 포함한다.
|
||||
|
||||
성공은 기존 **raw E residual, translation/rotation correction, chart/geometry, finite quantities, history/recovery gates**로 판정한다. `psi`만 작거나 transformed residual만 작다고 accept하지 않는다. 동일 trial에서 `v,a,A`를 재계산하고 기존 accepted frame identity에 맞춰 commit한다. 실패 trial은 폐기하며 history와 final output을 보존한다. Rechart·clamp·cutback·mass shift는 없다.
|
||||
|
||||
## 5. 일반 mesh와 부분 BC 경계
|
||||
|
||||
Element assembly 후에도 같은 nodal director/null convention으로 `N^T I=0`를 확인해야 한다. Mesh 연결로 새 unconstrained inertia null이 생기거나 geometry/rank 조건이 깨지면 현 후보의 nodal-rank 가정이 성립하지 않는다.
|
||||
|
||||
Zero fixed BC에서는 full variation을 `delta q=L delta q_f`로 놓는다. 필요한 basis는 `ker(M_ff)`의 **실제 admissible null**이며 nodal null에서 constrained 성분을 단순 제거한 벡터가 아니다. Current nodal block mass에서는 node의 constrained rotational component 중 null 성분이 하나라도 nonzero이면 그 nodal null은 free space에 속하지 않는다. Translation BC만으로는 rotational null이 제거되지 않는다. 이 구분은 chart에 따라 달라질 수 있어 support interval의 constant rank를 확인해야 한다.
|
||||
|
||||
General case는 `E_f=L^T E`, `N_f^T E_f=psi_f`, `T_f=[P_f;N_f^T]`가 가역인 local coordinate selection을 정의해야 한다. `Dpsi_f`에는 free-coordinate null derivative가 들어간다. 동적 null basis 선정/rank threshold/pivot 전략은 이번 all-free 고정-P probe로 승인되지 않는다.
|
||||
|
||||
현재 `phi_z` 제외 방식은 flat fixture의 `|n_z|`가 0에서 멀다는 조건에 한정한다. 일반 유한회전 mesh에서는 가장 큰 null 성분에 대응하는 coordinate row를 빼는 선택이 가능하지만, 결정적 pivot과 전환 시 local smoothness·conditioning을 독립 검토해야 한다. Silent fallback으로 `phi_z`를 항상 drilling coordinate라고 가정하지 않는다. Time-varying BC는 현재 제외된 범위이다.
|
||||
|
||||
## 6. 다른 후보와 조사 순서
|
||||
|
||||
| 후보 | 바뀌는 내용 | 조사상 의미와 한계 |
|
||||
| --- | --- | --- |
|
||||
| 본 `PE/psi` 정확 Newton | off-root algebraic equations의 linearization | Moving-null connection으로 인한 raw row 약화를 직접 분리. Endpoint physics/time scheme는 보존하나 nonlinear iteration 정책 채택 필요 |
|
||||
| Predictor 시작 | 첫 trial만 `q_n` 대신 extrapolated q 사용 | Off-root inertia를 줄일 수 있으나 원본 h=.1의 invalid geometry가 ADR-010의 채택 이유였다. 모든 지원 입력의 유효성 확인 없이는 전역 변경 불가 |
|
||||
| Line search/trust region | Newton step 크기 선택·trial sequence | Chart 이탈을 피할 수 있으나 이미 conditioning이 나쁜 방향의 구조를 직접 없애지 않는다. 단위에 맞는 merit function, trial rejection 및 기존 정책 변경 결정 필요 |
|
||||
| Physical/null condensation | reduced unknowns 또는 nested algebraic solve | `psi_q N` regularity와 implicit-function branch를 검증해야 한다. Moving basis derivatives·BC 처리가 필요하며 단순 raw matrix Schur solve는 원 correction을 바꾸지 않음 |
|
||||
| Rate/acceleration constraint recovery | discrete history 또는 별도 출력 | C08a/C08b 시간정확도 문제와 연결되지만 현재 runtime Newton 후보와 별도 문제. ADR-011의 미채택 recovery를 자동 승계하지 않음 |
|
||||
|
||||
## 7. 필수 반증과 미해결 결정
|
||||
|
||||
1. 동일 실패 직전 state에서 `N^T I=0`, `N^T K_eff=Dpsi-DN^T E`, `Dpsi` directional FD와 null derivative를 독립 검증한다. 물리/director-null·numerical drill contribution을 구별한다.
|
||||
2. `T`의 가역성, raw/변환 residual root 동등성, scaled conditioning과 linear backward error를 대조한다. Near-root 두 Newton correction의 차이가 residual과 함께 줄어드는지 확인한다.
|
||||
3. 동일 start/h/coefficients/iteration cap/strict gates에서 fixture의 실패 증분과 전체 5 grids를 비교한다. Physics/reference tolerance나 accepted chart를 완화하지 않는다.
|
||||
4. C08a 및 C08b의 q/v/physical-a 시간정확도·hidden1/hidden2·energy/momentum을 별도로 판정한다. Runtime 완료만으로 이 항목들을 PASS로 바꾸지 않는다.
|
||||
5. Production 채택에는 independent review, all-free 이외 실제 BC/mesh 지원 범위 검증, original reference·affected regression, ADR-010과의 algorithm 계약 결정이 남는다. 본 문서는 그 승인이나 수치 통과를 하지 않는다.
|
||||
|
||||
본 owner는 자체 수식 유도와 기존 계약 대조만 수행했다. Build/test/solver 실행은 하지 않았고, 기존 JSON의 수치를 새 통과 evidence로 재분류하지 않았다.
|
||||
@@ -0,0 +1,200 @@
|
||||
# 대수좌표 제거 DAE 연구: 수식과 실행 결과
|
||||
|
||||
상태: 연구 후보·독립 검토 대기. Owner: formulation, 2026-10-12. [연구의 권장 후보](../../../docs/hht-implicit-dynamics/research.md#temporal-remediation-research)를 소형 문제로 비교했다. Production·tests·PRD·ADR·reference·tolerance는 수정하지 않았다. 실제 shell 검증이나 새 적분법 채택이 아니다.
|
||||
|
||||
재현 자료는 [Decimal80 probe](compare.py), [실제 native 기록](run-native.json), [stdout](run.log), [전체 grid·오차·조건 결과](comparison.json)다. `python -B compare.py` 실행 native0이며 FESA build/CTest는 실행하지 않았다. 모든 문제는 `t∈[0,1]`, `h=1/16,...,1/256`, exact frame0의 q/v/a와 `A0=a0`, `alpha=-.05,beta=.275625,gamma=.55`를 사용했다. 초기 구간을 제외하지 않고 모든 `n+1` endpoint의 component infinity-norm 최대 오차를 구했다. 이 raw absolute error는 FESA N7의 normalized/noise-certified metric이 아니다.
|
||||
|
||||
## 1. 관찰과 범위
|
||||
|
||||
| 소형 문제·적분 | h=1/256 전체 최대 q / v / a 오차 | 마지막 q / v / a 차수 |
|
||||
| --- | --- | --- |
|
||||
| Cubic, full-coordinate | 0 / 1.532799e-5 / .032218337 | 미정 / 2 / 1 |
|
||||
| Cubic, reduced+미분 복구 | 0 / 0 / 0 | 오차0이므로 미정 |
|
||||
| Explicit-time, full-coordinate | 0 / 1.532799e-5 / .032218337 | 미정 / 2 / 1 |
|
||||
| Explicit-time, reduced+미분 복구 | 0 / 0 / 0 | 오차0이므로 미정 |
|
||||
| Moving-null oscillator, full-coordinate | 1.286321e-5 / 2.430511e-5 / .032224263 | 2.001052 / 2.000946 / 1.000776 |
|
||||
| Moving-null oscillator, reduced+미분 복구 | 1.977078e-6 / 3.548473e-6 / 9.930913e-6 | 1.999855 / 1.999676 / 1.999113 |
|
||||
| Moving-null+spectator BC | 각 적분의 바로 위 moving-null 결과와 동일 | 동일 |
|
||||
|
||||
Moving-null reduced 결과의 finest hidden0/1/2 최대는 각각 `2e-80/2e-80/2e-79`; 원 free endpoint E는 `1.4e-79` 이하다. Full-coordinate의 원 free E도 `3.22e-75` 이하지만 hidden1/2는 약 `1.53276e-5/.03222480`다. 따라서 이 probe에서는 endpoint 평형의 높은 정확도와 full-coordinate 가속도의 1차 startup 오차가 공존한다.
|
||||
|
||||
이 결과는 고정 null의 exact 해에 국한되지 않는 후보 근거다. 그러나 shell의 구성·하중·회전 branch·coupled modes·stiff frequency를 포함하지 않는다. 실제 C08a/C08b PASS, 일반 DAE 수렴 정리, production 알고리즘 승인을 주장하지 않는다. Decimal80 precision과 작은 root 잔차가 자체 독립 numerical review를 대신하지 않는다.
|
||||
|
||||
## 2. 공통 원 방정식과 국소 제거 계약
|
||||
|
||||
소형 모델은 dimensionless 좌표·시간·에너지를 사용하며 물리적 shell 단위나 재료모형을 대체하지 않는다. 원 covector residual은 `E(q,v,a,t)=I(q,v,a)+g(q,t)`, `I=M(q)a+c(q,v)`, `g=f_int-F`이다. 실제 admissible mass-null N에 대해 `N^T I=0`이고 원 대수 조건은
|
||||
|
||||
\[
|
||||
\psi(q,t)=N(q)^Tg(q,t)=0.
|
||||
\]
|
||||
|
||||
고정 BC를 정확히 lift한 free chart에서 `q=Q(p,t)`가 같은 ψ=0 branch를 parameterize한다고 가정한다. 독립 좌표 p의 dimension은 physical mass rank이다. `Q_p`가 full column rank이고 `[Q_p,N]`가 free space를 span해야 한다. 이때 ψ=0와
|
||||
|
||||
\[
|
||||
\boxed{R(p,\dot p,\ddot p,t)=Q_p^TE(Q,Q_p\dot p+Q_t,a,t)=0}
|
||||
\]
|
||||
|
||||
은 원 free E=0와 동등하다. 새로운 holonomic reaction multiplier를 넣거나 ψ를 다른 constraint로 교체하지 않는다. `M_r=Q_p^TMQ_p`는 positive definite이어야 한다. Numerical rank·regularity threshold는 이 문서에서 만들지 않는다.
|
||||
|
||||
독립 좌표에서만 기존 HHT의 physical/auxiliary acceleration recurrence를 적용한다.
|
||||
|
||||
\[
|
||||
\mathcal A_{n+1}=(1+\alpha)\ddot p_{n+1}-\alpha\ddot p_n,
|
||||
\quad p_{n+1}=p_n+h\dot p_n+h^2[(1/2-\beta)\mathcal A_n+\beta\mathcal A_{n+1}],
|
||||
\]
|
||||
\[
|
||||
\dot p_{n+1}=\dot p_n+h[(1-\gamma)\mathcal A_n+\gamma\mathcal A_{n+1}].
|
||||
\]
|
||||
|
||||
Full 좌표와 물리적 rate는 완전한 chain rule로 복구한다.
|
||||
|
||||
\[
|
||||
\boxed{q=Q(p,t),\quad v=Q_p\dot p+Q_t,\quad
|
||||
a=Q_p\ddot p+Q_{pp}[\dot p,\dot p]+2Q_{pt}\dot p+Q_{tt}.}
|
||||
\]
|
||||
|
||||
따라서 매끄러운 정확 embedding에서는 hidden0/1/2가 해석적으로 0이다. Full q/v/a는 일반적으로 기존 component별 Newmark recurrence를 만족하지 않는다. 이는 output-only null-a recovery와 다른 새 시간적분 계약이다. Full auxiliary A/history에 affine mapping으로 reduced `mathcal A`를 대입해서는 안 된다. 공개 state/schema와 production history 선택은 미결정이다.
|
||||
|
||||
국소 split `(p,z)`에는 같은 branch의 `G=psi_z` 가역성이 필요하다. `z_p=-G^-1 psi_p`, `z_t=-G^-1 psi_t`이며 고차 미분은 ψ=0를 미분해 구한다. `psi_q N` 가역성은 transversality 근거이지만 선택한 전역 `(p,z)` chart를 제공하지 않는다. Basis·branch·patch 전환에는 계약이 필요하며 constant rank, bounded derivative, M_r positivity가 전 구간에서 유지되어야 한다.
|
||||
|
||||
## 3. 정확한 reduced tangent 인계
|
||||
|
||||
고정 endpoint t에서 `c_v=gamma/(beta h)`, `c_a=1/[(1+alpha)beta h²]`다. HHT 종속 reduced rate와 Newton 좌표 변화 ξ에 대해
|
||||
|
||||
\[
|
||||
\delta q=Q_p\xi,\quad
|
||||
\delta v=Q_{pp}[\xi,\dot p]+Q_pc_v\xi+Q_{pt}\xi,
|
||||
\]
|
||||
\[
|
||||
\begin{aligned}
|
||||
\delta a={}&Q_{pp}[\xi,\ddot p]+Q_pc_a\xi+Q_{ppp}[\xi,\dot p,\dot p]
|
||||
+2Q_{pp}[c_v\xi,\dot p]\\
|
||||
&+2Q_{ppt}[\xi,\dot p]+2Q_{pt}c_v\xi+Q_{ptt}\xi.
|
||||
\end{aligned}
|
||||
\]
|
||||
\[
|
||||
\boxed{\delta R=(DQ_p[\xi])^TE+
|
||||
Q_p^T\{(D_qI+D_qg)\delta q+D_vI\delta v+M\delta a\}.}
|
||||
\]
|
||||
|
||||
`DQ_p` 항은 평형 밖에서 0이 아니다. 기존 내력 tangent가 낮은 차수의 미분을 쓴다고 embedding의 3차 미분과 explicit-time mixed 항을 생략할 수 없다. Shell 구현 전 analytic implicit differentiation 또는 독립 검증한 미분 경로가 필요하다. Exact frame0 q/v/a는 같은 branch와 미분 일관성 조건에 있어야 하며 reduced 초기값으로 그 값을 정확히 복구해야 한다.
|
||||
|
||||
## 4. 실험 모델과 full-coordinate 비교의 의미
|
||||
|
||||
**Cubic (NR-009 재현).** `q=(x,z)`, `T=xdot²/2`, `U=(z-x³)²/2`, `F=0`, `M=diag(1,0)`, `N=e_z`다. Exact 해는 `x=t,z=t³`, 초기값은 `x0=z0=0,vx0=1,vz0=0,a0=0`이다. Embedding `Q(p)=(p,p³)`에서 `M_r=1`, physical reduced equation은 `pddot=0`이다. Full-coordinate 비교는 exact endpoint q가 E=0을 만족하되 x와 z 모두 기존 Newmark rate를 유지한다. 기존 가속도 최대가 increment10에서 나타나는 결과를 재현했으며 finest에서는 t=10/256이다.
|
||||
|
||||
**Explicit time.** `T=xdot²/2`, synthetic `U(q,t)=(z-tx)²/2`, `F=(2,0)`, exact `x=t²,z=t³`, `a0=(2,0)`다. 시간 의존 potential 연구 모델이며 FESA 현 force-only 공개 입력의 확장이 아니다. `Q(p,t)=(p,tp)`이므로 `vz=p+t pdot`, `az=2 pdot+t pddot`다. Explicit/mixed-time 항을 빼면 oracle이 무효다. `M_r=1`, `pddot=2`이다.
|
||||
|
||||
**Moving-null oscillator.** `y=x+z²/2`, `T=ydot²/2`로 놓으면
|
||||
|
||||
\[
|
||||
U=\tfrac12y^2+\tfrac12(z-y^3)^2,\quad
|
||||
H=(1,z),\quad M=H^TH,\quad N=(-z,1)^T,
|
||||
\quad\psi=N^T\nabla U=z-y^3.
|
||||
\]
|
||||
\[
|
||||
I=H^T(a_x+za_z+v_z^2),\quad
|
||||
E=H^T(\ddot y+y)+(z-y^3)\nabla(z-y^3).
|
||||
\]
|
||||
|
||||
Exact reduced 해는 `y=sin(t),ydot=cos(t),yddot=-sin(t)`이며 Decimal80 sin/cos는 수렴하는 Taylor series로 계산했다. `Q(y)=(y-y⁶/2,y³)`는 `H Q_y=1`이므로 해석적으로 `M_r=1`, `MN=0`이다. Null 방향이 이동하며 그 미분과 embedding 곡률은 0이 아니다. 계산에서는 실제 full M과 N을 곱해 MN을 평가하고 `Q_y^T M Q_y`를 계산한다. M_r의 80자리 반올림 편차는 raw JSON에 보존한다.
|
||||
|
||||
**Production Newton의 수렴과 시간적분을 분리해 비교하기 위해** full-coordinate moving-null 쪽은 정확 endpoint-root oracle을 사용했다. Root finder에서 ψ=0을 대수적으로 대입하고 `a_x(q)+z a_z(q)+v_z(q)²+y=0`를 풀지만 **모든 full 좌표의 v/a/A는 기존 recurrence와 history**로 계산한다. 위치만 대수 제거한 root solve는 chain-rule rate를 쓰지 않으므로 reduced time integrator와 다르다. 이 residual/root 결과는 production Newton trial 경로를 검증하지 않는다.
|
||||
|
||||
**Spectator BC의 제한된 범위.** 독립 w에 `T_w=wdot²` (mass2), force3와 `w=0` BC를 추가했다. 최종 probe는 실제 3성분 q/v/a/A, full 3×3 M, full E를 구성하며 `L=[e_x,e_z]`로 free equation을 선택하고 constrained residual에서 reaction을 구한다. `E_w=R_w=-3`이고 원 full E max3은 그대로 보존하며 free E와 `E-R`은 near-zero다. Admissible null은 `(-z,1,0)`, reduced mass는1이다. 이는 분리된 spectator의 projection·reaction 검사일 뿐 **moving-null을 가르는 shell partial rotational BC의 검증이 아니다.** Rank 전환·구속 성분 변경도 검사하지 않았다. 최초 하드코딩 bookkeeping 버전은 `v1-bookkeeping-*` 원본 기록으로 보존하고 현재 검증과 구분한다.
|
||||
|
||||
## 5. 다음 검증과 미결정
|
||||
|
||||
- 독립 reviewer가 Decimal80 recurrence, moving-null full E, embedding 미분, exact startup과 all-time 차수를 재현해야 한다. 저장된 native0은 실행 완료만 뜻한다.
|
||||
- 실제 shell에는 locally integrable physical chart, 대수 branch solve, 하중의 explicit-time 의존성을 포함한 Q 미분, positive reduced mass와 admissible partial-BC null이 필요하다. 소형 모델 결과로 충분하지 않다.
|
||||
- 같은 계수와 exact frame0을 쓴다고 reduced history가 승인된 full-coordinate history와 같아지지 않는다. Production 채택에는 수식·수치 검토와 명시적인 algorithm/state 결정이 필요하다. Tolerance/reference 변경은 제안하지 않는다.
|
||||
- Full E·반력·hidden0/1/2·물리적 a와 공간 rate·energy/P/L·현행 C08·원본 비교는 각각 필수 evidence다. 이 oscillator는 frequency1이므로 C08b의 stiff-mode 시간 해상도 문제를 해결하지 않는다.
|
||||
|
||||
## 6. 실제 shell chart 후보와 구현-ready capability gap
|
||||
|
||||
현재 source는 읽기만 했다. [shell evaluator API](../../../include/fesa/elements/nonlinear_mitc4_shell.h)와 [내부 Jet·scalar constraint](../../../src/fesa/elements/nonlinear_mitc4_shell.cpp), [관성 API](../../../include/fesa/elements/nonlinear_mitc4_inertia.h)를 대조했다. 소형 prototype가 구현한 것은 해석적으로 알려진 polynomial Q와 독립 scalar reduced HHT뿐이다. Actual shell algebraic branch solver·Q 미분·reduced trial/state·assembly·tangent는 **구현하지 않았다**.
|
||||
|
||||
### 6.1 Chart 선택과 물리 동등성의 조건
|
||||
|
||||
| 후보 | 필요한 국소 조건 | 관성과 실제 구현 영향 |
|
||||
| --- | --- | --- |
|
||||
| `p=(u3,phi_x,phi_y)`, 대수 `z=phi_z` | Assembled coupled `G=psi_z` 가역성, 같은 branch, `psi_q N` 가역성, positive M_r 및 기존 chart/geometry | p의 회전2개는 director의 physical coordinates가 아니다. Algebraic branch에 따라 d와 M_r가 달라질 수 있으나 전체 Q/곡률/convective terms를 유지하면 원 full E와 동등할 수 있음 |
|
||||
| Initial-director 기준 η2와 fiber spin σ | Smooth director patch, swing·principal Log 및 spin branch 유효성, `psi_sigma` 가역성, positive physical mass | Director mapping은 σ와 무관하여 kinetic fiber를 명시적으로 분리. Rotation-vector mapping·미분·patch 전환의 새 계약 필요 |
|
||||
|
||||
첫 후보는 모든 node의 φz를 **전역 coupled vector z**로 대수 제거한다. G를 각 node의 scalar derivative로 대신하지 않는다. `G=psi_z` 가역성만으로 `[Q_p,N]` 가역성이 따라오지는 않는다. `ker(M)=span(N)`과 `psi_q Q_p=0` 아래에서는 M_r의 positivity가 `range(Q_p)∩span(N)={0}`를 확인하며 이 transversality가 원 잔차 동등성에 필요하다. `psi_q N` 가역성이 같은 조건을 표현한다.
|
||||
|
||||
C08b의 현재 수용된 회전이 약 .27rad라는 사실은 locality 근거이지 위 G/positivity의 증명이 아니다. Actual trajectory의 처음·최대·끝·원 실패 부근에서 root branch, G, M_r, `[Q_p,N]`를 확인하면 그 fixture의 bounded local prototype를 시작할 근거가 될 수 있다. 전 구간·일반 mesh·부분 회전 BC·다른 하중까지 승인되는 것은 아니다.
|
||||
|
||||
둘째 후보의 예시는 `d=eta_1 e_1+eta_2 e_2+sqrt(1-|eta|²)d0`, `d·d0>0`의 초기 director hemisphere와
|
||||
|
||||
\[
|
||||
R(d,\sigma)=R_{swing}(d;d_0)\exp([d_0]_\times\sigma),\qquad
|
||||
\phi=\operatorname{Log}R
|
||||
\]
|
||||
|
||||
이다. `R_swing d0=d`이므로 d·시간 rate와 kinetic mapping은 σ와 무관하다. 이 frame의 local M_r는 대수 spin branch와 무관하지만 η에 따른 configuration dependence와 관성 곡률은 남는다. Existing principal rotation/swing/twist 범위를 통과해야 하며 η hemisphere 또는 Log branch가 바뀌면 silent rechart할 수 없다. 고정 generalized force를 이 chart로 옮길 때에는 원 covector를 `Q_p^T F`로 pullback한다. 임의의 spatial torque로 재해석하지 않는다.
|
||||
|
||||
### 6.2 실제 보유 operator와 부족한 operator
|
||||
|
||||
| 수식 operator | 현재 확인된 capability | Reduced exact full Newton의 gap |
|
||||
| --- | --- | --- |
|
||||
| g, Dg | Internal force와 K_int, 외력은 별도 owner | 기존 force/assembly 의미 보존하여 사용 가능 |
|
||||
| D²g[v,v] | Optional `coordinate_direction`의 `force_directional_second`; U의3차 contraction | Mixed D²g[p,r]는 polarization으로 구할 수 있지만 비용·산술을 검토해야 함. D³g는 API에 없음 |
|
||||
| Internal ψ, Dψ, D²ψ | `drilling_constraint`, gradient, Hessian에 **raw γ가 아닌** `k_de*S_i*gamma_ei`를 저장 | Shared node assembly·동일 unit n 확인 필요. Full ψ에는 external null work가 추가됨 |
|
||||
| D³ψ 방향 contractions | Scalar Jet 내부 `mixed_second[j]`는 `D³scalar[e_j,v,v]`를 수송하지만 constraint용 반환 field에는 노출되지 않음 | 실제 scalar ψ에 대한 third contraction과 mixed polarization의 의미·finite 여부·독립 검산이 필요. Production API 추가는 미승인 |
|
||||
| n, Dn, D²n[v,v] | `DynamicRotationKinematics`의 unit null, gradient, directional_second | Force가 있는 ψ의 D³n와 일반 mixed second contraction 공개 operator가 없음 |
|
||||
| I, DqI, DvI, M | 기존 nonlinear inertia evaluation | Q·rate를 적용한 전체 chain rule와 DQ_p 항을 새 reduced assembly에서 연결해야 함 |
|
||||
| Force의 시간미분 | 현재 C08a oracle는 Fdot/Fddot | Mixed-time third ψ에는 F'''까지 필요할 수 있음. 현 oracle를 F'''가 이미 있는 것으로 보고하지 않음 |
|
||||
|
||||
일반적인 `psi=N^Tg`의 직접 미분은
|
||||
|
||||
\[
|
||||
D\psi[p]=(DN[p])^Tg+N^TDg[p],
|
||||
\]
|
||||
\[
|
||||
D^2\psi[p,r]=(D^2N[p,r])^Tg+(DN[p])^TDg[r]
|
||||
+(DN[r])^TDg[p]+N^TD^2g[p,r],
|
||||
\]
|
||||
\[
|
||||
\begin{aligned}
|
||||
D^3\psi[p,r,s]={}&(D^3N[p,r,s])^Tg
|
||||
+\sum_{sym}(D^2N[p,r])^TDg[s]\\
|
||||
&+\sum_{sym}(DN[p])^TD^2g[r,s]+N^TD^3g[p,r,s].
|
||||
\end{aligned}
|
||||
\]
|
||||
|
||||
마지막 식의 각 `sum_sym`은 서로 다른3항의 합이다. 이 경로는 D³g, 즉 내부 energy의4차 미분까지 요구한다. 그러나 **현 shell의 정확 scalar identity로 더 낮은 차수 경로가 가능하다.** Source는 `n_i^T grad(U_phys)=0`, `Dgamma_ej[n_i]=delta_ij*S_i`, `S_i=1/|J_l^-1 d_i|`를 사용하여
|
||||
|
||||
\[
|
||||
\psi_i(q,t)=\sum_{e\ni i}k_{de}S_i(q)\gamma_{ei}(q)-n_i(q)^TF_i(t)
|
||||
\]
|
||||
|
||||
를 직접 평가한다. 이 identity와 assembled nodal null convention을 독립 검증하면 scalar ψ의3차 미분으로 충분하며 generic U의4차 tensor를 만들 필요가 없다. C08b의 F=0에서는 ψ의 내부 scalar 3차 contraction이 핵심 gap이다. API의 `drilling_constraint_gradient/Hessian`를 raw γ의 미분이라 오해하지 않는다. 원 ψ의 단순화는 내력에서 U_drill를 제거하는 근거가 아니다.
|
||||
|
||||
F(t)가 q-independent이고 n은 q-only인 C08a convention에서는
|
||||
|
||||
\[
|
||||
\psi_{qt}[p]=-Dn[p]^TF',\quad
|
||||
\psi_{qqt}[p,r]=-D^2n[p,r]^TF',\quad
|
||||
\psi_{qtt}[p]=-Dn[p]^TF'',\quad \psi_{ttt}=-n^TF'''.
|
||||
\]
|
||||
|
||||
Pure q의 third 항에는 `-D³n^T F`도 필요하다. Constant public CLOAD에서는 F'/F''/F'''가0이지만 F가 nonzero이면 pure-q null 미분항은 사라지지 않는다. Synthetic time-dependent potential에는 위 force-only 단축식 대신 전체 ψ 미분을 쓴다.
|
||||
|
||||
### 6.3 Implicit embedding 미분의 최소 인계
|
||||
|
||||
`tilde Q(p,t)=(Q(p,t),t)`로 time을 추가한 chart를 쓰면 pure/mixed derivatives를 같은 규칙으로 처리할 수 있다. `A=D tilde Q[a]`, `B=D tilde Q[b]`, `C=D tilde Q[c]`라 놓으면 대수 좌표의2차·3차 방향 미분은
|
||||
|
||||
\[
|
||||
G\,D^2z[a,b]=-D^2\psi[A,B],
|
||||
\]
|
||||
\[
|
||||
G\,D^3z[a,b,c]=-D^3\psi[A,B,C]
|
||||
-\sum_{sym}D^2\psi[D^2\tilde Q[a,b],C].
|
||||
\]
|
||||
|
||||
p/t의 affine identity 부분은2차 이상에서0이며 RHS는 이전 단계의 Q 미분을 포함한다. 이 규칙은 Q_pp, Q_pt, Q_tt와 δR에 필요한 Q_ppp/Q_ppt/Q_ptt를 제공한다. Dense global3·4차 tensor 전체를 요구하는 것은 아니지만, 실제 필요한 mixed directional contractions와 같은 G의 정확한 solve가 필요하다. Perturbed algebraic root solve를 사용한 FD reduced tangent는 **private diagnostic** 경로로 연구할 수 있다. 이는 production analytic δR의 capability gap을 닫거나 새 수치기준을 승인하지 않는다.
|
||||
|
||||
부분 BC에서는 먼저 exact zero-constrained lift를 정의하고 실제 admissible null의 dimension을 정해야 한다. Constraint가 nodal null을 가르면 단순히 그 null의 constrained 성분을 삭제해서는 안 된다. p/z 수와 G 크기는 실제 rank와 함께 정하며 rank·pivot·chart 변경 규약을 정하지 않은 채 all-free5/node를 보편적으로 쓰지 않는다. Empty-null은 기존 full-coordinate 문제로 환원된다.
|
||||
|
||||
이 후보가 C08b의 큰 physical-complement fixed-grid 오차를 고친다는 보장은 없다. Actual prototype에는 current5grids의 full q/v/physical-a 및 공간 rates, 초기/최대/최종 오차 위치, dynamic spectrum·mode amplitude·omega*h, 원 free E/reaction, hidden0/1/2, energy/P/L을 함께 비교하는 조건을 남긴다. 국소 제거의 null-rate 개선과 고주파 temporal resolution의 개선 여부를 별도로 판단해야 한다.
|
||||
@@ -0,0 +1,64 @@
|
||||
# ADR-018 graph BL/solve certificate controller draft
|
||||
|
||||
Status: private design only. Shared test source is frozen while the configured
|
||||
Debug focused run is active.
|
||||
|
||||
## Scope
|
||||
|
||||
Stored `ConstraintRowOperators::{gradient,hessian,time_gradient,third_directional_gradient}`
|
||||
are exact binary64 inputs. The draft certifies only fixed-basis contraction and
|
||||
linear graph-solve arithmetic. It does not certify AD/FE operator formation,
|
||||
root branch sensitivity, or nonlinear history.
|
||||
|
||||
## Test-local data
|
||||
|
||||
`GraphArithmeticCertificate` should retain dense-small matrices/vectors:
|
||||
`Ghat`, `P_hat`, `BG`, `BP`, `bhat`, `Bb`, `R`, `Hbar`, `Q`, `rbar`, `Bx`.
|
||||
It is private to `hht_shell_physics_test.cpp`.
|
||||
|
||||
For every graph row and basis column, replay source order:
|
||||
|
||||
1. Iterate `row.gradient`, then `column.support`; match full DOF.
|
||||
2. Product radius is `ProductRadius(C, 0, L, BL)`.
|
||||
3. Reduce matching products in `RateSumRadius` source order.
|
||||
|
||||
This yields stored `Ghat=C*Lz`, `P_hat=C*Lp`, and entrywise `BG/BP`.
|
||||
The same replay builds RHS radii for `P*vp`, `P*ap`, Hessian/time contractions.
|
||||
No zero-BL fallback is allowed.
|
||||
|
||||
## Inverse action
|
||||
|
||||
For each canonical unit RHS, use the existing graph solver and N2-L check to
|
||||
obtain one stored column of `R`. Replay `I-Ghat*R` with product/sum radii to
|
||||
make nonnegative `Hbar`. Require finite values and `thetaH=||Hbar||inf<1`.
|
||||
Construct
|
||||
|
||||
`Q[i,j] = abs(R[i,j]) + (abs(R)*Hbar*1)[i] / lower(1-thetaH)`.
|
||||
|
||||
All bound additions/divisions use existing outward helpers. Failure of the
|
||||
positive lower denominator or either theta is INCONCLUSIVE, not a radius reset.
|
||||
|
||||
## Uncertain solve
|
||||
|
||||
Replay `rhat=bhat-Ghat*xhat`; combine absolute residual and rounding into
|
||||
`rbar`. With actual RHS bound `Bb`:
|
||||
|
||||
`y = Q*(rbar + Bb + BG*abs(xhat))`, `A=Q*BG`, `thetaG=||A||inf`.
|
||||
|
||||
Require `thetaG<1`; set
|
||||
|
||||
`Bx[i] = y[i] + rowSum(A,i)*normInf(y)/lower(1-thetaG)`.
|
||||
|
||||
Undo any actual rate scaling before using `Bx` in physical lifted radii.
|
||||
|
||||
## Tests
|
||||
|
||||
Add only after the current focused run closes:
|
||||
- small chart with plane coefficient value zero / positive BL, proving BG/BP
|
||||
remain nonzero after an omitted-BL mutation;
|
||||
- 1x1 and 4x4 source-order graph replay, N2-L checked;
|
||||
- one thetaH and one thetaG failure must return incomplete;
|
||||
- use a known 4x4 perturbation from the independent 32-case probe as a
|
||||
behavior oracle, without copying implementation equations as its oracle.
|
||||
|
||||
No production API, tolerance, grid, or sparse global inverse is added.
|
||||
@@ -0,0 +1,215 @@
|
||||
# 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).
|
||||
- 관성의 전개식·초기 일관성·정확도 근거는 [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#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)는 완료되었으며 표시된 수식에서 새 오류는 발견되지 않았다. 질량 선택은 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-probes.md)은 variable-mass SDOF와 moving-null 예제에서 추천 후보의 약 2차 q/v/physical-a 수렴을 관찰했고 naive 확대의 약 1차 오차를 구별했다. 이는 shell solver 검증이나 일반 drilling DAE 증명을 대신하지 않는다. Build/test·solver/reference 실행은 문서 작업 범위가 아니므로 미실행이다. 독립 검산과 Numerical Review 결과는 해당 owner의 보고서가 소유하며 본 문서가 스스로 승인하지 않는다.
|
||||
@@ -0,0 +1,209 @@
|
||||
# 유한회전 관성·초기화·HHT 상세
|
||||
|
||||
[주 정식화](formulation.md)의 구현 수식 상세다. 상태와 미결정 항목은 주 문서를 따른다. 아래 전개는 명시한 FESA 운동학에서의 자체 유도이며 Abaqus 내부 절차의 기술이 아니다.
|
||||
|
||||
사용자 선택에 따라 이번 HHT 기본 및 구현 질량은 consistent `T_C/M_C`다. Lumped 식은 비교용 이론·기존 eigen 회귀 참고로 보존하며 이번 HHT 구현·검증 요구를 추가하지 않는다. 다른 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는 시간에 따라 갱신하지 않는다.
|
||||
|
||||
비교용 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 보정항을 추가하지 않는다.
|
||||
|
||||
## 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은 비교식에만 사용한다. 긴 유도는 선행 [mass 계약 M3](../mitc4-shell-feast-eigen/formulation-mass.md#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 일반식으로 사용하지 않는다. 사용자는 보간된 shell 속도의 운동에너지를 직접 적분하는 consistent 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-probes.md)의 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을 제공하지 않으므로 현재 구현이 이 알고리즘을 지원한다고 보고하지 않는다.
|
||||
@@ -0,0 +1,79 @@
|
||||
# C08b linear momentum: 정확한 이산 불변량과 판정 후보
|
||||
|
||||
상태: 수식 검토 후보·독립 numerical review 및 사용자 결정 대기. Owner: formulation, 2026-10-12. 공식 acceptance·tolerance·tests·production은 변경하지 않았다. 현재 C08b FAIL을 유지한다. 이 문서는 새 수치 tolerance를 정하지 않는다.
|
||||
|
||||
기준은 [현재 시간적분 정식화](../../../docs/hht-implicit-dynamics/formulation.md), [lumped 관성 상세](../../../docs/hht-implicit-dynamics/formulation-inertia.md), [현재 연구](../../../docs/hht-implicit-dynamics/research.md#temporal-remediation-research)다. 별도 [DAE 비교](formulation.md)는 full-coordinate와 reduced 시간적분의 구별을 소유한다. 아래 증명은 **현재 C08b의 원 full-coordinate HHT**에 대한 자체 유도이며 reduced 후보를 자동 포함하지 않는다.
|
||||
|
||||
## 1. 가정과 검산
|
||||
|
||||
적용 대상은 flat, centered, all-free C08b; constant reference nodal kinetic measure; 모든 node의 scalar first thickness moment `s_i=0`; 외력·반력 없음이다. 이 경우 병진 질량은 고정 `m_i I3`이고 병진 관성 covector는 정확히 `I_ui=m_i a_ui`다. 회전 관성·convective 항은 병진 행에 결합하지 않는다. 물리적 linear momentum은
|
||||
|
||||
\[
|
||||
P_n=\sum_i m_i v_{ui,n},\qquad [P]=mass\cdot length/time.
|
||||
\]
|
||||
|
||||
Objective physical/drilling potential은 균일 translation `u_i<-u_i+c`에 불변이다. `c`에 대한 directional derivative를 취하면
|
||||
|
||||
\[
|
||||
0=D U[\delta u_i=c]=c\cdot\sum_i f_{int,ui}
|
||||
\quad\Longrightarrow\quad\sum_i f_{int,ui}=0.
|
||||
\]
|
||||
|
||||
따라서 매 endpoint의 원 방정식 `E_ui=m_i a_ui+f_int,ui=0`를 합하면
|
||||
|
||||
\[
|
||||
S_{a,n}:=\sum_i m_i a_{ui,n}=0.
|
||||
\]
|
||||
|
||||
초기화가 원 dynamic equilibrium을 만족하면 `S_a,0=0`이며 `A0=a0`에서 `S_A,0=0`이다. 고정된 m_i를 HHT 필터 식과 velocity update에 그대로 적용하면
|
||||
|
||||
\[
|
||||
S_{A,n+1}:=\sum_i m_i A_{ui,n+1}
|
||||
=(1+\alpha)S_{a,n+1}-\alpha S_{a,n}=0,
|
||||
\]
|
||||
\[
|
||||
P_{n+1}-P_n=h\{(1-\gamma)S_{A,n}+\gamma S_{A,n+1}\}=0.
|
||||
\]
|
||||
|
||||
결론은 **유효한 endpoint 해를 수용한 모든 h에서 exact arithmetic의 P_n=P_0**이다. P_0가0일 필요가 없고 fixed-step의 h 크기·alpha 부호도 이 algebraic 보존 증명에 영향을 주지 않는다. 이 증명은 그 h의 비선형 해가 존재하거나 Newton이 수렴한다는 보장은 아니다. Total translation displacement도 `Σm_i u_i,n=Σm_i u_i,0+t_n P_0`로 affine하게 진행한다.
|
||||
|
||||
## 2. Strict temporal decrease와 정확 branch의 충돌
|
||||
|
||||
각 grid의 error를 기존 방식으로 `e_P(h)=max_n |P_n-P_0|` 또는 같은 norm/scale의 비음수 metric으로 평가하면 exact branch에서 모든 h의 e_P는0이다. 따라서
|
||||
|
||||
\[
|
||||
e_P(h/2)<e_P(h)\quad\Longrightarrow\quad0<0
|
||||
\]
|
||||
|
||||
는 거짓이다. Ratio `e_P(h/2)/e_P(h)`도0/0으로 정의되지 않는다. P를 정확히 보존하는 구현이 strict-decrease를 반드시 통과할 수 있다는 계약은 이 조건 아래 모순이다. 이 결론은 실제 double drift가0이라는 뜻이나 현재 threshold를 면제한다는 뜻은 아니다.
|
||||
|
||||
현재 binary의 비영 drift는 finite residual·내력 translation-identity 산술오차·q−q_p 차감·rate/history update·momentum 합산의 noise 또는 구현 오류일 수 있다. 이 각각의 원인을 구분해야 하며 보존 정리만으로 현재 관찰값을 roundoff라고 인증하지 않는다. 다만 원 full-coordinate HHT의 정확 endpoint branch에는 **P의 시간절단 오차가 없으므로**, 작은 nonzero drift의 refinement ratio를 HHT temporal order라고 해석할 수 없다.
|
||||
|
||||
## 3. 현재 bound를 보존하는 진단 인계
|
||||
|
||||
새 threshold 없이 다음 raw quantity를 같은 state와 같은 nodal measure로 평가할 수 있다.
|
||||
|
||||
\[
|
||||
S_a=\sum_i m_i a_{ui},\quad S_A=\sum_i m_i A_{ui},\quad
|
||||
F_{sum}=\sum_i f_{int,ui},\quad E_{sum}=\sum_i E_{ui},
|
||||
\]
|
||||
\[
|
||||
r_A=S_{A,n+1}-(1+\alpha)S_{a,n+1}+\alpha S_{a,n},
|
||||
\]
|
||||
\[
|
||||
r_v=P_{n+1}-P_n-h[(1-\gamma)S_{A,n}+\gamma S_{A,n+1}].
|
||||
\]
|
||||
|
||||
원 balance에서는 `S_a=E_sum-F_sum`이다. `r_A/r_v`는 실제 계산한 full rates/history가 이산 관계와 맞는지 검사하며, finite-precision residual/identity/update/aggregation 결함을 분리하는 근거다. Independent high-precision·compensated sum 대조와 실제 native/source identity는 실행 담당자의 별도 evidence로 남긴다. 이 문서 owner는 추가 binary/probe 실행을 하지 않았다.
|
||||
|
||||
**계약 검토 후보:** 승인된 현재 absolute P bound `<=1e-3`를 유지하고 P에는 절대 보존·원 balance·update consistency를 검사한다. Temporal-decrease는 이 특수 문제에서 P를 제외하고 기존 E/L quantity에 유지하는 선택을 검토한다. E/L의 기존 acceptance를 새로 증명하거나 완화하는 제안은 아니다. Precision/noise 문제의 처리와 bound 적용 의미는 numerical/I/O review 후 사용자에게 구체적으로 제시해야 한다.
|
||||
|
||||
사용자는 연구의 권장 작업 순서를 승인했으며 **이 acceptance 변경을 승인하지 않았다.** 따라서 지금 tests/공식 checklist를 수정하거나 P FAIL을 PASS로 재분류해서는 안 된다.
|
||||
|
||||
## 4. 적용 한계
|
||||
|
||||
- `s_i!=0`인 일반 geometry에서는 `P=Σ(m_i v_ui+s_i B_i v_phi,i)`이고 구성·속도에 따른 항이 존재한다. 위 constant translational block 증명을 그대로 쓸 수 없다.
|
||||
- 외력·지지 반력이 있으면 P 변화는 net impulse와 대조해야 하며0 보존 검사를 적용하지 않는다. Constrained equation을 free sum에 넣어 보존을 만들지 않는다.
|
||||
- Translation invariance는 U_drill까지 포함한 실제 discretization의 identity여야 한다. 회전 객관성만으로 구현된 translation sum의 정확성을 대체하지 않는다.
|
||||
- Angular momentum L에는 current lever arm·director rate가 들어가며 위 선형 velocity 합산 증명과 다르다. Continuous Noether balance가 HHT의 exact discrete L 보존을 보장하지 않는다. HHT energy도 별도 numerical dissipation·work 검사 대상이다.
|
||||
- Reduced-coordinate 후보에서 recovered full translation rate가 같은 원 Newmark recurrence를 따르지 않는다면 이 증명을 자동 승계할 수 없다. 실제 reduced residual/translation symmetry·update의 독립 불변량 검산이 필요하다.
|
||||
Reference in New Issue
Block a user