feat(linear-static-3d-euler-beam): step 3 - numerical-review

2026-08-09 02:52:47 +09:00
parent 06b285a8ba
commit 2bccffb8c5
1 changed files with 228 additions and 0 deletions
@@ -0,0 +1,228 @@
 # Linear Static 3D Euler Beam Numerical Review
 ## Metadata
 - feature_id: `linear-static-3d-euler-beam`
 - source_formulation: `docs/formulations/3d-isoparametric-euler-beam-formulation.md`
 - source_requirements: `docs/requirements/linear-static-3d-euler-beam.md`
 - source_research: `docs/research/linear-static-3d-euler-beam-research.md`
 - approved_design: `docs/superpowers/specs/2026-08-08-linear-static-3d-euler-beam-design.md`
 - status: `pass-for-implementation-planning`
 - owner_agent: `numerical-review-agent`
 - date: `2026-08-09`
 - scope: formulation correctness and implementation-planning readiness only
 ## Review Verdict
 - verdict: `pass-for-implementation-planning`
 - reason: The DOF/sign convention, natural boundary actions, 4x12 B matrix, 12x12 stiffness, 2-point Gauss rule, transformation, constant line-load vector, constrained partition, and residual reaction are mutually consistent. No confirmed mathematical defect was found.
 - qualification: This verdict permits implementation planning only; it is not reference-comparison, physics-sanity, release, or full Abaqus-compatibility approval.
 ## Critical Findings
 ### Confirmed defects
 - None.
 ### Independent DOF, curvature, and end-action derivation
 The fixed local nodal order is
 $$
 [u_1,v_1,w_1,\theta_{x1},\theta_{y1},\theta_{z1},
 u_2,v_2,w_2,\theta_{x2},\theta_{y2},\theta_{z2}].
 $$
 For a small right-handed section rotation,
 $\boldsymbol\theta\times[0,y,z]^T
 =[z\theta_y-y\theta_z,-z\theta_x,y\theta_x]^T$.
 The zero-shear constraints are therefore
 $\theta_z=v'$ and $\theta_y=-w'$. Hence
 $$
 \varepsilon_{xx}=u'-zw''-yv'
 =\varepsilon_0+z\kappa_y-y\kappa_z,
 \quad \kappa_y=-w'',\quad \kappa_z=v''.
 $$
 Using the reviewed section-force definitions independently gives
 $$
 M_y=EI_y\kappa_y=-EI_yw'',\qquad
 M_z=EI_z\kappa_z=EI_zv''.
 $$
 Twice integrating the bending virtual work, with outward sign $n=-1$ at
 $x=0$ and $n=+1$ at $x=L$, yields
 $$
 F_y=-nM_z',\quad M_z^{end}=nM_z,\qquad
 F_z=nM_y',\quad M_y^{end}=nM_y.
 $$
 The axial and torsional pairs are $F_x=nN$ and $M_x=nT$. Direct expansion of
 $K_ld_l$ reproduces these relations at both endpoints when the V0 distributed
 load is zero. Thus the equilibrium end action is an outward action, while
 $DBd$ is the positive-local-x section resultant before the outward-normal
 sign is applied. These outputs must remain distinct.
 ### B matrix, stiffness, quadrature, and modes
 - Differentiating the Hermite fields reproduces the documented 4x12 B matrix. The $\kappa_y=-w''$ row has $[-H_1'',H_2'',-H_3'',H_4'']$ in the $[w_1,\theta_{y1},w_2,\theta_{y2}]$ columns; the $\kappa_z=v''$ row has $[H_1'',H_2'',H_3'',H_4'']$ in the $[v_1,\theta_{z1},v_2,\theta_{z2}]$ columns.
 - Axial/torsional B rows are constant and bending rows are linear in $\xi$, so $B^TDBJ_x$ is degree at most two. The 2-point Gauss rule, exact through degree three, exactly produces the closed-form 12x12 stiffness.
 - The axial and torsion blocks each have rank one, and each exact bending block has rank two. With positive $EA,GJ,EI_y,EI_z$, the free-element stiffness is symmetric positive semidefinite with rank 6.
 - The six rigid modes are the three equal nodal translations, equal $\theta_x$, $\theta_{y1}=\theta_{y2}$ with $w_2-w_1=-L\theta_y$, and $\theta_{z1}=\theta_{z2}$ with $v_2-v_1=L\theta_z$. Each gives $Bd=0$.
 - Since $d^TKd=\int(Bd)^TD(Bd)dx$, energy is strictly positive outside the six-dimensional rigid kernel, but the unconstrained element itself is not positive definite.
 - A 1-point bending rule reduces the total element rank to 4 and introduces two spurious zero-energy modes. Production must use the approved 2-point rule or the reviewed closed form.
 ### Transformation and energy invariance
 Projection/normalization of `n1` gives orthogonal unit $e_x,e_y$, and
 $e_z=e_x\times e_y$ makes $R=[e_x^T;e_y^T;e_z^T]$ orthogonal and
 right-handed: $RR^T=I$, $\det R=+1$. With
 $T=\operatorname{diag}(R,R,R,R)$ in node translation/rotation block order,
 $$
 d_l=Td_g,\quad K_g=T^TK_lT,\quad f_g=T^Tf_l,
 \quad d_g^TK_gd_g=d_l^TK_ld_l.
 $$
 This also confirms virtual-work invariance. Reversing the transform direction
 would violate the documented component contract.
 ### Constant local line-load vector
 Exact integration of $N_q^T[p_x,p_y,p_z,m_x,0,0]^T$ gives equal half-span
 nodal forces/torques and
 $$
 M_{y1}=-p_zL^2/12,\quad M_{y2}=+p_zL^2/12,\qquad
 M_{z1}=+p_yL^2/12,\quad M_{z2}=-p_yL^2/12.
 $$
 The signs follow from the $-H_2,-H_4$ interpolation for $w$ and
 $H_2,H_4$ for $v$. The 2-point rule integrates the cubic load integrand
 exactly. This remains a formulation-only kernel: `*DLOAD` input, a Domain
 distributed-load object, and CLI distributed loading are out of scope.
 ### Independent numerical cross-check
 An independent double-precision evaluation at $L=3.7$ with unequal positive
 $EA,GJ,EI_y,EI_z$ produced:
 | check | observed result |
 | --- | ---: |
 | 2-point Gauss/closed-form normalized error | $3.20\times10^{-16}$ |
 | symmetry normalized error | $2.00\times10^{-17}$ |
 | six rigid-mode residual norms | 0 to $2.81\times10^{-15}$ |
 | numerical rank at relative $10^{-10}$ | rank 6 |
 | constant line-load normalized error | $3.08\times10^{-16}$ |
 | 1-point stiffness | rank 4 |
 | rotated local/global energy difference | 0 |
 | $\det R$ | $0.9999999999999999$ |
 These calculations corroborate the derivation but do not replace production
 tests.
 ## Numerical Risk Assessment
 Confirmed defects, risks, and open issues are separated.
 | id | risk or limit | assessment and required control |
 | --- | --- | --- |
 | NR-R01 | zero/near-zero length | Singular Jacobian and divergent $L^{-3}$ terms; `NR-T05` tests the approved scale-aware boundary before evaluation. |
 | NR-R02 | zero/tangent-parallel `n1` | Undefined axes and inertia directions; `NR-T05` tests the projection boundary and forbids a silent fallback. |
 | NR-R03 | under-integration | 1-point integration adds two mechanisms; `NR-T03` must show 2-point rank 6 and negative-control rank 4. |
 | NR-R04 | mixed translation/rotation scaling | Raw rank is length-unit dependent; `NR-T04` uses length scaling and a well-scaled fixture. |
 | NR-R05 | incomplete constraints | `Kff` is singular while any rigid mode survives; `NR-T09` covers free, partial, and stable systems. |
 | NR-R06 | extreme stiffness ratios/units | A mathematically SPD system may be ill-conditioned; `NR-T10` checks residual/analytical error. Automatic scaling or a new threshold is not approved. |
 | NR-R07 | transform direction | Rotated bending signs/components can be wrong; `NR-T05` checks orthogonality, handedness, work, and energy. |
 | NR-R08 | end-action/section-cut confusion | Opposite signs can contaminate comparison; `NR-T07` checks both endpoints and distinct result types. |
 | NR-R09 | sparse duplicate reduction | Thread-dependent summation can break reproducibility; `NR-T11` checks canonical CSR and normalized values. |
 | NR-R10 | short/deep beam | Euler–Bernoulli may be too stiff; explicit limitation only. No unapproved slenderness diagnostic and no B31 reinterpretation. |
 | NR-R11 | warping, eccentricity, or `I12!=0` | The diagonal constitutive model is insufficient; explicit rejection/exclusion, never silent approximation. |
 | NR-R12 | result averaging | Averaging can hide discontinuities; baseline nodal averaging is excluded and `NR-T07` preserves endpoint identity. |
 - shear_locking: N/A because no independent shear strain exists; deep-beam applicability remains limited.
 - volumetric_locking: N/A.
 - hourglass: N/A for the approved 2-point rule; the distinct 1-point rank deficiency is NR-R03.
 - distortion: no curved/distorted mapping exists in this straight affine element; length and axis singularities are NR-R01/NR-R02.
 - nonlinear convergence: N/A; linear residual and conditioning are NR-T09/NR-T10.
 - mesh convergence: tip-force/tip-moment polynomial fields are one-element exact; the formulation-only constant transverse load has a quartic exact displacement and requires refinement evidence.
 ## Consistency Checks
 | check | result | evidence |
 | --- | --- | --- |
 | units | pass | Stiffness blocks and nodal work pairs are dimensionally consistent. |
 | dimensions | pass | B is 4x12, D is 4x4, and $B^TDBJ_x$ is 12x12. |
 | signs | pass | Cross-product, curvature, natural-boundary, end-action, and line-load derivations agree. |
 | dof_ordering | pass | Node blocks preserve `[u,v,w,theta_x,theta_y,theta_z]`. |
 | coordinate_transforms | pass | R/T are right-handed orthogonal maps with local = transform times global. |
 | jacobian | pass with validation | $J_x=L/2>0$ for accepted elements. |
 | constitutive | pass | Positive diagonal $EA,GJ,EI_y,EI_z$ gives nonnegative energy; `I12!=0` is rejected. |
 | integration | pass | The 2-point locations, unit weights, and $J_x$ integrate stiffness and constant-load kernels exactly. |
 | symmetry | pass | Variational and numerical checks agree. |
 | rank 6 and rigid modes | pass | Independent ranks $1+1+2+2=6$ and six rigid vectors span the kernel. |
 | positive deformation energy | pass | Strictly positive outside the rigid kernel. |
 | constrained partition | pass | $K_{ff}d_f=F_f-K_{fc}d_c$ includes nonzero $d_c$. |
 | residual reaction | pass | $R_c=(Kd-F)_c=K_{cf}d_f+K_{cc}d_c-F_c$. |
 | output locations | pass with downstream contract | Gauss values, section resultants, end actions, and residual reactions remain distinct. |
 ## Verification Readiness
 | test id | concrete criterion |
 | --- | --- |
 | `NR-T01-DOF-SIGN` | Compare `B*d` with independently differentiated Hermite fields at multiple $\xi$; normalized error $\le10^{-12}$ and explicit $\theta_y=-w'$, $M_y=-EI_yw''$. |
 | `NR-T02-PATCH` | Check endpoint value/slope, constant axial strain/twist, both constant-curvature patches, and six rigid modes; rigid residual $\le10^{-10}$. |
 | `NR-T03-GAUSS-CLOSED` | All 12x12 entries and symmetry normalized $\le10^{-12}$; negative 1-point control must be rank 4. |
 | `NR-T04-RANK-ENERGY` | Use $Q=\operatorname{diag}(1,1,1,L,L,L,1,1,1,L,L,L)$ and $\widehat K=Q^{-T}KQ^{-1}$. For a documented well-scaled fixture require six singular values $\le10^{-10}\sigma_{max}$, six $>10^{-10}\sigma_{max}$, and positive energy for six deformation vectors. |
 | `NR-T05-TRANSFORM` | For a non-axis-aligned member require normalized orthogonality, handedness, virtual-work, transformed-stiffness, and energy errors $\le10^{-12}$; test geometry thresholds. |
 | `NR-T06-LINE-LOAD` | Compare all 12 signed closed-form components at normalized $\le10^{-12}$; separately verify `*DLOAD` remains unsupported. |
 | `NR-T07-END-SIGNS` | Pure axial/torsion/two-plane bending states must match outward natural actions and positive-face resultants at both endpoints, normalized $\le10^{-12}$; no averaging. |
 | `NR-T08-ANALYTICAL` | Axial, torsion, local-y, and local-z cantilevers at relative $\le10^{-9}$; constant transverse load requires load-vector equality and mesh convergence. |
 | `NR-T09-PARTITION-SPD` | Nonzero $d_c,F_c$: displacement and free/reaction residual normalized $\le10^{-10}$. Stable `Kff` factorizes; surviving rigid modes fail structurally with no result. |
 | `NR-T10-CONDITIONING` | Sweep documented valid scales; accepted well-conditioned cases require residual $\le10^{-10}$ and analytical error $\le10^{-9}$; numerically unresolved cases fail explicitly. |
 | `NR-T11-DETERMINISTIC-ASSEMBLY` | Across repeated thread counts require identical CSR structure/reduction order and values normalized $\le10^{-12}$. |
 ### Open issues
 1. `NR-O01-PARDISO-CONTRACT`: official oneMKL evidence is still required for SPD matrix type, phases, zero-based CSR, repeated RHS, conditioning/error reporting, and failure codes before that adapter plan is approved. This is non-blocking for element mathematics and connects to NR-T09/NR-T10.
 2. `NR-O02-DETERMINISTIC-REDUCTION`: stable COO sort and duplicate-summation rules are project policy and must be made explicit before NR-T11.
 3. `NR-O03-STATION-NORMALIZATION`: reversed connectivity/local-axis orientation and legitimate jumps at loaded interior nodes need an explicit downstream row-normalization/eligibility rule. The legacy baseline may use its documented stable orientation and unloaded interior stations, but mismatch must never be averaged. NR-T07 covers element signs.
 No open issue requires formulation revision. NR-O01/NR-O02 are implementation-planning
 handoffs; NR-O03 belongs to I/O and reference-model contracts.
 ## Required Revisions
 ### Formulation Agent
 - None. No confirmed mathematical defect or missing derivation blocks implementation planning.
 ### Research Agent
 - Resolve NR-O01 from official Intel documentation before finalizing the PARDISO adapter plan.
 ### Reference Model Agent
 - Make NR-O03 orientation and unloaded-interior assumptions explicit without modifying the approved legacy artifacts.
 ## Downstream Handoff
 ### Implementation Planning Agent
 - Convert NR-T01 through NR-T11 into TDD `RED -> GREEN -> VERIFY` work without changing approved tolerances.
 - Preserve stiffness/partition/factorization before load/effective-RHS/substitution.
 - Treat SPD as a post-constraint property and retain singularity/conditioning diagnostics.
 - Resolve NR-O01 and NR-O02 before the corresponding implementation steps.
 ### I/O Definition Agent
 - Preserve outward endpoint action, positive-face section resultant, Gauss result, and assembled residual reaction as distinct identities.
 - Define NR-O03 without broadening the approved Abaqus subset.
 ### Reference Model Agent
 - Cover axial, torsion, both bending planes, rotated-space transformation, endpoint signs, and nonzero prescribed displacement.
 - Keep `reference/cantilever beam/` read-only, stress comparison N/A, and `*DLOAD` outside CLI support.