From 2bccffb8c5aa4a350c0c81d70f0b1b23738ed390 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Sun, 9 Aug 2026 02:52:47 +0900 Subject: [PATCH] feat(linear-static-3d-euler-beam): step 3 - numerical-review --- .../linear-static-3d-euler-beam-review.md | 228 ++++++++++++++++++ 1 file changed, 228 insertions(+) create mode 100644 docs/numerical-reviews/linear-static-3d-euler-beam-review.md diff --git a/docs/numerical-reviews/linear-static-3d-euler-beam-review.md b/docs/numerical-reviews/linear-static-3d-euler-beam-review.md new file mode 100644 index 0000000..74a3e35 --- /dev/null +++ b/docs/numerical-reviews/linear-static-3d-euler-beam-review.md @@ -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.