diff --git a/docs/formulations/mitc4-shell-formulation.md b/docs/formulations/mitc4-shell-formulation.md index 3520521..93d6851 100644 --- a/docs/formulations/mitc4-shell-formulation.md +++ b/docs/formulations/mitc4-shell-formulation.md @@ -5,9 +5,13 @@ - feature_id: `linear-static-mitc4-shell` - source_requirement: `docs/requirements/linear-static-mitc4-shell.md` - source_research: `docs/research/linear-static-mitc4-shell-research.md` +- source_numerical_review: `docs/numerical-reviews/linear-static-mitc4-shell-review.md` - status: `ready-for-numerical-review` - owner_agent: `formulation-agent` - date: `2026-08-11` +- revision_basis: numerical review commit `0a5aad4`; findings `NR-C01` through + `NR-C05` and decisions `NR-D01` through `NR-D02` +- revision_state: `ready-for-numerical-rereview-not-implementation-planning` - product_execution_scope: `small-strain, small-rotation linear static only` - future_formulation_scope: `geometrically nonlinear Total Lagrangian residual/tangent; not executable` @@ -37,10 +41,11 @@ Source labels `S4` and `S4R` both select this one FESA formulation by an approve semantic mapping. They do not select Abaqus integration or stabilization behavior, and this document makes no Abaqus formulation-equivalence claim. -Exact drilling scale/coefficient, drilling-energy warning ratio, smooth-director -angle, drilling-load projection tolerance, and Jacobian/warpage thresholds remain -Numerical Review decisions. Symbols are retained for them; no implementation default -is implied. +Exact drilling reference family/coefficient, drilling-energy warning ratio, +smooth-director angle, and geometry thresholds remain Numerical Review decisions. +The common dimensionless drilling coordinate, geometry-measure inventory, +drilling-load projection tolerance, and normalized algebraic checks are fixed below +by the first Numerical Review. No remaining open symbol is an implementation default. ## 2. Scope and assumptions @@ -192,7 +197,20 @@ using the same four midsurface Gauss locations as the element formulation. ### 4.2 Smooth shared-node director -For source node `I` and its incident accepted shell elements `E(I)`, define +For source node `I`, sort its incident accepted shell elements `E(I)` by stable +source-element identity. Before any averaging, require every unordered incident +pair to have the same orientation hemisphere: + +$$ +c_{pair,I}=\min_{e0. +$$ + +For a node with only one incident element this pair test is vacuous. A nonpositive +pair dot product is an `opposed-incident-normal` geometry error; it is not repaired +by flipping a candidate or by letting cancellation occur in the average. This +pairwise rule is deterministic and stronger than testing only the final average. + +After the pairwise test, define $$ \mathbf s_I=\sum_{e\in E(I)} A_e\mathbf n_e, @@ -200,8 +218,8 @@ $$ \mathbf d_I=\frac{\mathbf s_I}{\|\mathbf s_I\|}. $$ -Incident elements are accumulated in stable source-element order. Before averaging, -opposing orientations fail. After averaging, every incident deviation +Incident elements are accumulated in the same stable order. After averaging, every +incident deviation $$ \theta_{eI}=\cos^{-1}\!\left( @@ -352,12 +370,19 @@ Prescribed global shell DOFs form the Dirichlet data. The only external generali forces are aggregated nodal `CLOAD` values. For a nodal moment `M_I`, define $$ -M_{d,I}=\mathbf d_I\cdot\mathbf M_I. +\rho_{M,I}=\frac{|\mathbf d_I\cdot\mathbf M_I|}{\|\mathbf M_I\|}. $$ -If `M_d,I` is nonzero beyond an approved scale-aware projection tolerance, model -validation fails as `unsupported-drilling-load`. The tolerance is a Numerical -Review/I-O decision and is not replaced by an arbitrary absolute epsilon. +An exactly zero aggregated moment is accepted as a separate zero case and does not +form this ratio. Every nonzero moment is accepted only when + +$$ +\boxed{\rho_{M,I}\le 10^{-12}}. +$$ + +Otherwise model validation fails as `unsupported-drilling-load`. The denominator is +never clamped. This is a direction/projection check and is invariant under a +consistent moment-unit conversion. Distributed natural-boundary terms may appear in the parent continuum derivation, but their existence is not an executable load contract. @@ -366,31 +391,43 @@ but their existence is not an executable load contract. ### 7.1 Physical virtual work -Let `V_20` be the admissible five-DOF element trial space and `V_20^0` its homogeneous -test space. The physical linear equilibrium is +Let `V_24` be the admissible global six-DOF element trial space after essential +boundary data are applied, and let `V_24^0` be its homogeneous test space. With $$ -\delta W_{int}^{phys}-\delta W_{ext}=0 -\qquad\forall\delta\mathbf q_{20}\in V_{20}^0, +\mathbf q_{20}=\mathbf T_p\mathbf q_g, +\qquad +\delta\mathbf q_{20}=\mathbf T_p\delta\mathbf q_g, +$$ + +the MITC-projected physical internal virtual work is + +$$ +\delta W_{int}^{phys} +=\int_{\Omega_0} +(\delta\overline{\boldsymbol\epsilon})^T +\mathbf C_5 +\overline{\boldsymbol\epsilon}\,dV +=\delta\mathbf q_{20}^T\mathbf f_{int}^{20} +=\delta\mathbf q_g^T\mathbf T_p^T\mathbf f_{int}^{20}, $$ where $$ -\delta W_{int}^{phys} -=\int_{\Omega_0} -\delta\overline{\boldsymbol\epsilon}^{,T} -\mathbf C_5 -\overline{\boldsymbol\epsilon}\,dV, +\mathbf f_{int}^{20} +=\int_{\Omega_0}\overline{\mathbf B}^T +\mathbf C_5\overline{\boldsymbol\epsilon}\,dV. $$ -and `overbar` denotes the MITC-projected shear field. In the current product path +The overbar denotes the common MITC projection used by strain, residual, and +stiffness. The current product external virtual work is $$ \delta W_{ext}=\delta\mathbf q_g^T\mathbf f_e^{CLOAD}. $$ -### 7.2 Numerical drilling potential +### 7.2 Numerical drilling potential and complete weak form Drilling regularization is algebraically separate: @@ -404,6 +441,22 @@ $$ It does not approximate a shell strain energy and produces no physical generalized strain, resultant, or stress. +The complete linear 24-DOF equilibrium statement is therefore + +$$ +\boxed{ +\delta\mathbf q_g^T\left[ +\mathbf T_p^T\mathbf f_{int}^{20} ++\mathbf T_d^T\mathbf K_d^l\boldsymbol\gamma +-\mathbf f_e^{CLOAD} +\right]=0 +\qquad\forall\delta\mathbf q_g\in V_{24}^0.} +$$ + +This is the normative weak form. A physical-only statement in `V_20` is merely the +restriction of this equation and must not be combined with a 24-DOF external-work +term as if both variations occupied the same test space. + ## 8. Bilinear discretization and degenerated geometry ### 8.1 Shape functions @@ -517,22 +570,88 @@ $$ - any additional bottom/middle/top recovery evaluation location before recovery is committed. +The center `(0,0,0)` is also included in the orientation and warpage inventory. Duplicate nodes, self-intersection, degenerate midsurface area, and reversed -connectivity are separate fail-closed geometry errors. +connectivity are separate fail-closed geometry errors. No failed location is +discarded or replaced by a value from another point. ### 9.3 Scale-aware measures -A dimensionless local determinant measure is +Let the consecutive midsurface edge inventory be + +$$ +\mathcal E=\{(1,2),(2,3),(3,4),(4,1)\}, +\qquad +\boxed{L_e=\max_{(I,J)\in\mathcal E}\|\mathbf X_J-\mathbf X_I\|}. +$$ + +`L_e` must be finite and strictly positive. It is the common element length used by +geometry checks and the DOF scaling in Section 12.5. At every distinct in-plane +location in the center, Gauss, tying, and committed recovery inventory, define + +$$ +a(\xi,\eta)=\|\mathbf A_\xi\times\mathbf A_\eta\|, +\qquad +\boxed{a_g(\xi,\eta)=\frac{a(\xi,\eta)}{L_e^2}}, +$$ + +$$ +\mathbf n_s(\xi,\eta)= +\frac{\mathbf A_\xi\times\mathbf A_\eta}{a(\xi,\eta)}, +\qquad +c_d(\xi,\eta)=\mathbf n_s\cdot\overline{\mathbf d}. +$$ + +`a_g` is the normalized surface-collapse/aspect measure. In particular, for +`A_xi=(1,0,0)` and `A_eta=(0,epsilon,0)` with `L_e=O(1)`, `a_g -> 0` as +`epsilon -> 0`; the angular measure below alone cannot detect that collapse. + +At every full three-dimensional validation point, define the dimensionless +angular/director determinant measure $$ j_s=\frac{J} {\|\mathbf G_\xi\|\,\|\mathbf G_\eta\|\,\|\mathbf G_\zeta\|}. $$ -A valid point requires `J>0` and `j_s>tau_J`; a valid element also satisfies an -approved distortion/warpage inventory such as `J_min/J_max`. The exact `tau_J` and -warpage bounds are Numerical Review decisions. No `max(1, geometry_scale)` clamp is -permitted because it destroys consistent-unit scaling. +For all positive finite point determinants, define the element-variation measure + +$$ +\boxed{r_J=\frac{J_{min}}{J_{max}}}, +\qquad +J_{min}=\min_{p\in\mathcal P_V}J_p, +\quad +J_{max}=\max_{p\in\mathcal P_V}J_p, +$$ + +and the surface-normal warpage measure relative to the center normal + +$$ +\boxed{\theta_w= +\max_{p\in\mathcal P_S} +\cos^{-1}\!\left(\operatorname{clamp} +(\mathbf n_s(0,0)\cdot\mathbf n_s(p),-1,1)\right)}. +$$ + +Here `P_V` contains every volume Gauss point, tying point at `zeta=0`, center, and +every committed bottom/middle/top recovery point; `P_S` contains their distinct +in-plane projections. A valid element must satisfy, without denominator clamping, + +$$ +\boxed{ +J_p>0,\quad +j_{s,p}>\tau_{ang},\quad +a_{g,p}>\tau_{area},\quad +c_{d,p}>\tau_{dir},\quad +r_J>\tau_{var},\quad +\theta_w<\theta_{warp}.} +$$ + +The measures and their location inventory are fixed by this formulation revision. +The positive dimensionless thresholds remain `needs-numerical-calibration`; they +must separate valid distortion/warp sweeps from collapsed negative sequences before +Numerical Review may approve them. No `max(1, geometry_scale)`, zero denominator, +failed-point omission, or pointwise orientation repair is permitted. ## 10. Linear kinematics and MITC4 shear projection @@ -736,15 +855,36 @@ $$ \qquad k_{d,I}>0, $$ -with rotational-stiffness units `force*length`. Write +with rotational-stiffness units `force*length`. The common physical normalization is $$ -k_{d,I}=\alpha_d k_{ref,I}, +D_{iso}=\frac{Et^3}{12(1-\nu^2)}, +\qquad +\boxed{\rho_{d,I}=\frac{k_{d,I}}{D_{iso}}}. $$ -where both the dimensionless coefficient `alpha_d` and the reference family remain -unapproved until Numerical Review. Dimensionally compatible candidates carried from -research are: +`rho_d,I` is dimensionless and is the only common coordinate for comparing drilling +families. For any candidate `c` written as + +$$ +k_{d,I}^{(c)}=\alpha_d^{(c)}k_{ref,I}^{(c)}, +$$ + +the candidate-specific conversion is + +$$ +\boxed{ +\rho_{d,I}^{(c)}=\alpha_d^{(c)} +\frac{k_{ref,I}^{(c)}}{D_{iso}}, +\qquad +\alpha_{d,I}^{eq,(c)}=\rho_{d,I}^{(c)} +\frac{D_{iso}}{k_{ref,I}^{(c)}}.} +$$ + +The second expression is the nodewise equivalent coefficient for a target `rho_d,I`. +A single actual candidate coefficient may therefore generate a range of `rho_d,I`; that +entire range is part of the calibration evidence. The dimensionally compatible +candidate distributions carried from research are: 1. transverse-shear/area transition family @@ -765,15 +905,38 @@ $$ 2. isotropic bending rigidity $$ -k_{ref}^{(B)}=D_{iso}=\frac{Et^3}{12(1-\nu^2)}; +k_{ref}^{(B)}=D_{iso}, +\qquad +\rho_{d,I}^{(B)}=\alpha_d^{(B)}; $$ 3. a documented positive statistic formed only from the physical rotational block - of `K_20`, whose entries all have `force*length` units. + of `K_20`, whose entries all have `force*length` units, converted by the same + `k_ref/D_iso` ratio. A raw statistic is not comparable until this conversion is + reported. + +The first two candidate scales differ sharply in the thin-shell limit: + +$$ +\boxed{\displaystyle \lim_{A_{eI}/t^2\to\infty} +k_{ref,I}^{(A)}/D_{iso}=6(1-\nu)/q}. +$$ + +For `nu=0.3` and `q=2.5e-5`, this ratio is `168000`. Consequently the same raw +coefficient, including `10^-3`, cannot represent the same small drilling stiffness +for candidates A and B. The thesis rule `10^-3 min(all K_ii)` is not admissible because it can mix -translational `force/length` and rotational `force*length` diagonals. The order -`10^-3` may only center a sensitivity study; it is not a nominal coefficient. +translational `force/length` and rotational `force*length` diagonals. A sweep must +instead expand logarithmically in actual `rho_d,I` until it brackets both: + +1. a low-side scaled-rank/conditioning or factorization failure; and +2. a high-side physical `U/N/M/Q` contamination boundary. + +A nominal value may be proposed only as the smallest point in a stable plateau, with +the adjacent lower and higher decades and separate physical/drilling energies +reported. The reference family, plateau, nominal value, and response/energy bounds +remain `needs-numerical-calibration`; no common `10^-3` center is retained. ### 12.3 Stabilized 24-DOF matrix @@ -802,7 +965,9 @@ For a valid free isolated element, the expected physical rank is 14. Embedding i 24 coordinates creates the six physical rigid modes plus four drilling null modes. Four positive independent `k_d,I` values should remove only those drilling modes, giving expected stabilized rank 18 and nullity 6. These are verification targets, -not substitutes for a computed singular-value/rank study. +not substitutes for the scaled singular-value/rank study defined in Section 12.5. +The exact-arithmetic rank statement is independent of the calibrated numerical-rank +threshold. ### 12.4 Energy split @@ -827,6 +992,51 @@ The ratio `E_drill/E_phys` is reported only when mathematically classifiable. If `E_phys` is zero or near zero, the two energies are reported explicitly; no arbitrary denominator clamp is used. The warning ratio remains open. +### 12.5 DOF scaling for rank and conditioning evidence + +Raw shell stiffness matrices mix translational and rotational units and therefore +must not be used for singular-value, eigenvalue, rank, or condition-number +acceptance. Define diagnostic coordinate scalings + +$$ +\mathbf q_{20}=\mathbf S_{20}\widehat{\mathbf q}_{20},\qquad +\mathbf S_{20}=\operatorname{blockdiag}_{I=1}^{4}(L_e\mathbf I_3,\mathbf I_2), +$$ + +$$ +\mathbf q_g=\mathbf S_{24}\widehat{\mathbf q}_g,\qquad +\mathbf S_{24}=\operatorname{blockdiag}_{I=1}^{4}(L_e\mathbf I_3,\mathbf I_3). +$$ + +The corresponding scaled element matrices are + +$$ +\boxed{\widehat{\mathbf K}_{20}=\mathbf S_{20}^T\mathbf K_{20}\mathbf S_{20}, +\qquad \widehat{\mathbf K}_e=\mathbf S_{24}^T\mathbf K_e^{24}\mathbf S_{24}}. +$$ + +Every entry of these matrices has units `force*length`. Rigid and other diagnostic +vectors are transformed by `q_hat=S^{-1}q`. Scaling is used only for algebraic +evidence; it does not modify the physical assembly, prescribed values, or solve. + +For a global shell model, use the deterministic model length + +$$ +L_m=\max_{e\in\mathcal E_{active}}L_e +$$ + +and construct the analogous full-model scaling with `(L_m I_3,I_3)` per six-DOF +node. Restrict it to the stable free-DOF order as `S_f` and define + +$$ +\widehat{\mathbf K}_{ff}=\mathbf S_f^T\mathbf K_{ff}\mathbf S_f. +$$ + +Global condition and numerical-rank evidence uses `K_hat_ff`; a valid `0 x 0 Kff` +case is classified separately and is not reported as singular. Numerical +rank/condition thresholds remain calibration decisions, but no raw mixed-unit +matrix may be used to choose them. + ## 13. Numerical integration ### 13.1 Stiffness rule @@ -1019,15 +1229,17 @@ locations. This formulation does not define nodal averaging or extrapolation. ## 15. Future-only geometrically nonlinear residual and tangent > **Non-executable boundary:** This section records the requested future formulation. -> It must not be routed through the current linear-static procedure or state. +> It must not be routed through the current linear-static procedure or state. The +> equations through Section 15.3 live only in a 20-coordinate physical director +> chart; they are not a complete global 24-DOF nonlinear element. ### 15.1 Total Lagrangian kinematics -Let the physical nonlinear coordinate vector retain three translations and two -director parameters per node. In the reference configuration, +Let `q_20` be the physical nonlinear coordinate vector with three translations and +two director parameters per node. In the reference configuration, $$ -\mathbf x(\xi,\eta,\zeta;\mathbf q) +\mathbf x(\xi,\eta,\zeta;\mathbf q_{20}) =\sum_I N_I(\mathbf X_I+\mathbf u_I) +\frac{t\zeta}{2}\sum_I N_I\mathbf d_I(\boldsymbol\phi_I), $$ @@ -1050,13 +1262,13 @@ $$ \mathbf E=\frac12(\mathbf F^T\mathbf F-\mathbf I). $$ -For physical coordinates `q_A` and `q_B`, define +For physical chart coordinates `q_20,A` and `q_20,B`, define $$ -\mathbf h_A=\frac{\partial\mathbf x}{\partial q_A}, +\mathbf h_A=\frac{\partial\mathbf x}{\partial q_{20,A}}, \qquad \mathbf h_{AB}=\frac{\partial^2\mathbf x} -{\partial q_A\partial q_B}. +{\partial q_{20,A}\partial q_{20,B}}. $$ The covariant Green-Lagrange components and their exact first and second derivatives @@ -1124,8 +1336,8 @@ Use second Piola-Kirchhoff stress `S` work-conjugate to `E`. In engineering-vect notation, define $$ -\mathbf B_A(\mathbf q)= -\frac{\partial\overline{\mathbf e}}{\partial q_A}. +\mathbf B_A(\mathbf q_{20})= +\frac{\partial\overline{\mathbf e}}{\partial q_{20,A}}. $$ For the geometrically nonlinear but materially linear candidate, @@ -1145,17 +1357,14 @@ Then the physical internal residual is $$ \boxed{ -r_A^{phys}(\mathbf q) +r_A^{phys}(\mathbf q_{20}) =\int_{\Omega_0}\mathbf B_A^T\mathbf s\,dV}. $$ -The total residual is - -$$ -\mathbf r(\mathbf q)= -\mathbf r^{phys}(\mathbf q) -+\mathbf r^{drill}(\mathbf q)-\mathbf f_{ext}(\mathbf q). -$$ +Collecting these entries gives the physical chart residual +`r_20^phys(q_20) in R^20`. A total 24-DOF residual is deliberately not written at +this point: drilling and external load can be combined with the physical residual +only after the nonlinear 20-to-24 coordinate map in Section 15.4 is selected. ### 15.3 Consistent material and geometric tangent @@ -1166,10 +1375,10 @@ $$ {\partial\overline{\mathbf e}}, \qquad \mathbf G_{AB}=\frac{\partial^2\overline{\mathbf e}} -{\partial q_A\partial q_B}. +{\partial q_{20,A}\partial q_{20,B}}. $$ -The consistent physical tangent is +The consistent `20 x 20` physical-chart tangent is $$ \boxed{ @@ -1189,56 +1398,94 @@ K_{AB}^{geo}=\int_{\Omega_0} $$ `K_geo` contains the stress-dependent second variation of Green-Lagrange strain, -including consistent finite-director derivatives. The future total tangent is +including consistent finite-director derivatives. These equations define +`K_20^phys=K_20^mat+K_20^geo`; they do not yet define congruence to a global +24-DOF tangent. For a constant generalized dead-load vector in this selected +20-coordinate chart, its chart load tangent is zero. Follower pressure and other +configuration-dependent loads are outside scope. Finite-rotation nodal-moment work +must be separately defined before assuming its load tangent is zero. + +### 15.4 Unresolved 20-to-24 mapping and drilling boundary + +Let a future finite global six-DOF coordinate vector be `q_g` and suppose an +approved nonlinear map exists: $$ -\boxed{ -\mathbf K_T= -\mathbf K^{mat}+\mathbf K^{geo}+\mathbf K^{drill}-\mathbf K^{load}}. +\mathbf q_{20}=\boldsymbol\Phi(\mathbf q_g),\qquad +\mathbf A(\mathbf q_g)=\partial\boldsymbol\Phi/\partial\mathbf q_g. $$ -For a constant generalized dead-load vector in the selected coordinate chart, -`K_load=0`. Follower pressure and other configuration-dependent loads are outside -scope. Finite-rotation nodal-moment work must be separately defined before assuming -its load tangent is zero. - -### 15.4 Future drilling tangent boundary - -If the linear initial-frame penalty is formally retained, +Here `A` is a `20 x 24` mapping Jacobian. Conditional on that map, virtual work gives $$ -\Pi_d=\frac12\boldsymbol\gamma(\mathbf q)^T -\mathbf K_d^l(\mathbf q)\boldsymbol\gamma(\mathbf q), +\mathbf r_g^{phys}=\mathbf A^T\mathbf r_{20}^{phys}. +$$ + +Its consistent global physical tangent must contain both the congruence and the +coordinate-map curvature term: + +$$ +\boxed{\mathbf K_g^{phys}=\mathbf A^T\mathbf K_{20}^{phys}\mathbf A+\sum_{a=1}^{20}r_{20,a}^{phys}(\partial^2\Phi_a/\partial\mathbf q_g^2)}. +$$ + +The candidate two-parameter director chart in Section 15.1 does not by itself +define `Phi`: a finite three-component global rotation coordinate, its +director-parallel gauge, chart update/recentering, and the first and second +derivatives of the complete map are still unspecified. + +Likewise, a future drilling potential would have to be objective and expressed in +the same global coordinates: + +$$ +\Pi_d=\frac12\boldsymbol\gamma(\mathbf q_g)^T +\mathbf K_d^l(\mathbf q_g)\boldsymbol\gamma(\mathbf q_g), \qquad -\mathbf r^{drill}=\frac{\partial\Pi_d}{\partial\mathbf q}, +\mathbf r_g^{drill}=\frac{\partial\Pi_d}{\partial\mathbf q_g}, \qquad -\mathbf K^{drill}=\frac{\partial^2\Pi_d}{\partial\mathbf q^2}. +\mathbf K_g^{drill}=\frac{\partial^2\Pi_d}{\partial\mathbf q_g^2}. $$ -However, a constant initial-frame diagonal penalty is not automatically objective -under large rotation. The nonlinear frame update, drilling coordinate, scale update, -and their first/second derivatives are unresolved future decisions. They do not -block the current linear kernel but do block geometrically nonlinear implementation. +A constant initial-frame diagonal penalty is not automatically objective under +large rotation. The nonlinear frame, drilling coordinate, scale update, and all +their first/second derivatives remain unresolved. Until those choices and +configuration-dependent load work are approved, neither a complete global residual +nor a complete global 24-DOF tangent exists. This does not block the current linear +kernel but blocks geometrically nonlinear implementation. ### 15.5 Newton equation and linear limit -At iteration `k`, +Only after Section 15.4 is closed may the global residual and tangent be defined as $$ -\mathbf K_T(\mathbf q^{(k)})\Delta\mathbf q -=-\mathbf r(\mathbf q^{(k)}), -\qquad -\mathbf q^{(k+1)}=\mathbf q^{(k)}+\Delta\mathbf q. +\mathbf r_g=\mathbf A^T\mathbf r_{20}^{phys}+\mathbf r_g^{drill}-\mathbf f_{ext}, +\qquad \mathbf K_T=\mathbf K_g^{phys}+\mathbf K_g^{drill}-\mathbf K_g^{load}. $$ -At the undeformed, stress-free state with the small-rotation director derivative, +The conditional Newton equation at iteration `k` is $$ -\mathbf K^{geo}=\mathbf0, +\mathbf K_T(\mathbf q_g^{(k)})\Delta\mathbf q_g=-\mathbf r_g(\mathbf q_g^{(k)}), +\qquad \mathbf q_g^{(k+1)}=\mathbf q_g^{(k)}+\Delta\mathbf q_g. $$ -and the physical material tangent reduces to `K_phys^24` in Section 12. With the -approved linear drill block, the complete limit is `K_e^24`. +At the undeformed, stress-free state, recovering the current linear formulation +requires all of the following identities: + +$$ +\mathbf A_0=\mathbf T_p,\qquad \mathbf r_{20}^{phys}=\mathbf0, +\qquad \mathbf K_{20}^{geo}=\mathbf0. +$$ + +Then the coordinate-map curvature term vanishes and + +$$ +\mathbf K_g^{phys}=\mathbf T_p^T\mathbf K_{20}\mathbf T_p=\mathbf K_{phys}^{24}. +$$ + +If the future map also recovers `gamma=T_d q_g` and the approved constant linear +drill block in this limit, the complete tangent reduces to `K_e^24`. This is a +conditional consistency requirement, not evidence that a nonlinear `Phi` or +objective drilling potential has already been selected. ## 16. Algorithm pseudocode @@ -1247,12 +1494,14 @@ approved linear drill block, the complete limit is `K_e^24`. ```text for each shell element in stable source order: validate four distinct source nodes and source-order geometry + compute L_e from the four consecutive midsurface edges compute center positive normal candidate n_e compute 2x2 surface area A_e for each shell source node in stable source order: gather incident candidates in stable element order - reject opposed, degenerate, or too-sharp incident normals + reject any nonpositive pairwise incident-normal dot product before averaging + reject degenerate or too-sharp incident normals d_I = normalize(sum(A_e * n_e)) select least-aligned global axis with deterministic tie break construct right-handed (a_I, b_I, d_I) @@ -1265,13 +1514,15 @@ input: X_I, d_I, a_I, b_I, E, nu, t, q_g, nodal CLOAD share build T, T_p, T_d initialize K20[20,20] = 0 +evaluate the complete center/Gauss/tying/recovery geometry inventory +validate pointwise J, j_s, a_g, c_d and aggregate r_J, theta_w evaluate and validate four midsurface tying locations for each 2x2 midsurface Gauss location in fixed order: construct and validate local frame (e1,e2,e3) compute the four direct covariant tying shear B columns at zeta=0 for each thickness Gauss point in fixed order: build degenerated geometry, J, reciprocal bases, H derivatives - validate finite positive scale-aware Jacobian + reuse the validated point geometry and reciprocal bases form direct covariant B columns replace only xi-zeta and eta-zeta shear by MITC tying interpolation transform to local engineering B_bar[5,20] @@ -1279,9 +1530,10 @@ for each 2x2 midsurface Gauss location in fixed order: check K20 finite and symmetric within approved normalized tolerance Kphys24 = T_p^T * K20 * T_p -construct candidate-aware positive Kd_local without closing open coefficient choices +construct positive Kd_local and report every k_d,I through rho_d,I = k_d,I/D_iso Kdrill24 = T_d^T * Kd_local * T_d Ke24 = Kphys24 + Kdrill24 +form S20, S24, Khat20, and Khat_e for rank/conditioning evidence only fint24 = Ke24 * q_g residual24 = fint24 - f_CLOAD return matrices, residual, transforms, frames, and separate energy operators @@ -1291,9 +1543,11 @@ return matrices, residual, transforms, frames, and separate energy operators ```text assemble all Ke24 contributions with stable element-local COO ordering +form model-length DOF scaling and Khat_ff for global rank evidence only partition full K into Kff, Kfc, Kcf, Kcc in stable free/constrained order factorize Kff before load assembly assemble and deterministically aggregate nodal CLOAD +accept an exact-zero nodal moment separately; otherwise require rho_M <= 1e-12 solve Kff * df = Ff - Kfc * dc reconstruct full displacement d compute full residual r = K*d - F @@ -1317,11 +1571,13 @@ compute E_physical and E_drill separately; never clamp a near-zero denominator ### 16.5 Future nonlinear tangent check ```text -given a future approved finite-director map and nonlinear drill potential: - evaluate residual r(q) - evaluate analytic consistent K_T(q) +given a future approved 20-coordinate finite-director chart: + evaluate physical residual r20(q20) and consistent K20(q20) + verify its directional derivative +given an approved global Phi map, objective drill potential, and load work: + evaluate global residual rg(qg) and complete consistent K_T(qg) for several normalized perturbation directions p and decreasing h: - compare K_T*p with [r(q+h*p)-r(q-h*p)]/(2*h) + compare K_T*p with [rg(qg+h*p)-rg(qg-h*p)]/(2*h) require the error to decrease in the expected truncation range ``` @@ -1331,18 +1587,42 @@ given a future approved finite-director map and nonlinear drill potential: - Shape-function partition of unity, Kronecker delta, and derivative sums. - Nodal and integration frames orthonormal and right-handed. -- Finite positive Jacobian at every Gauss, tying, and committed recovery point. -- `K20`, `Kphys24`, `Kdrill24`, and `Ke24` symmetry. +- The complete `J/j_s/a_g/c_d/r_J/theta_w` inventory at center, Gauss, tying, + and committed recovery points. +- Scaled `K20`, `Kphys24`, `Kdrill24`, and `Ke24` symmetry and spectrum. - Transformation work/energy invariance. -- Physical rigid modes satisfy normalized `Kphys24 r=0`. +- Physical rigid modes satisfy normalized scaled stiffness action. - Stabilized free-element nullity is exactly six; accepted non-rigid physical modes have positive physical energy. - Pure drill vectors have zero physical energy and positive drilling energy. - Consistent force/length unit rescaling leaves dimensionless decisions unchanged. -Requirements currently prescribe normalized targets `1e-12` for frame, symmetry, -and transformation-energy checks and `1e-10` for rigid action, residual, and global -equilibrium unless Numerical Review documents replacements. +For any nonzero scaled stiffness under test, the approved normalized checks are +the following. Here `r_hat` denotes a constructed scaled rigid-mode vector, not the +assembled residual. + +$$ +e_{sym}=\|\widehat{\mathbf K}-\widehat{\mathbf K}^T\|_F/\|\widehat{\mathbf K}\|_F\le10^{-12}, +$$ + +$$ +e_{rigid}=\|\widehat{\mathbf K}\widehat{\mathbf r}\|_2/(\|\widehat{\mathbf K}\|_2\|\widehat{\mathbf r}\|_2)\le10^{-10}, +$$ + +$$ +e_{frame}=\|\mathbf R^T\mathbf R-\mathbf I\|_F\le10^{-12}. +$$ + +For a nonzero transformation-energy case, + +$$ +e_T=|E_g-E_l|/(|E_g|+|E_l|)\le10^{-12}. +$$ + +A zero matrix/vector denominator is invalid test construction and is not clamped to +pass. When both transformation energies are exactly zero, the rigid or pure-null +case is classified by its separate stiffness-action test. Linear-system residual +and global-equilibrium evidence retain the normalized `1e-10` requirement. ### 17.2 Patch and sign checks @@ -1369,16 +1649,20 @@ Independently verify: ### 17.4 Drilling sensitivity -For every candidate reference scale, sweep dimensionless coefficients around the -research order of magnitude and record: +For every candidate reference scale, convert candidate coefficients to the actual +nodewise `rho_d,I` inventory. Expand a logarithmic sweep until both the low-side +rank/conditioning failure and high-side physical-response contamination boundary +are observed, and record: -- free-element rank and constrained system conditioning; +- free-element scaled rank and scaled constrained-system conditioning; - global `U` and physical `N/M/Q` sensitivity; - `E_phys` and `E_drill` without denominator clamping; - invariance under consistent unit conversion, thickness ratios, and mesh sizes. -The study must select the reference family, nominal coefficient/range, and warning -criterion before Implementation Planning. +Candidate runs are compared by overlapping actual `rho_d,I` ranges, never by equal +raw `alpha_d`. The study must select the reference family, the smallest stable +plateau value, its adjacent-decade sensitivity, and the energy warning criterion +before Implementation Planning. ### 17.5 Reference-comparison boundary @@ -1393,8 +1677,10 @@ not expected to be numerically identical on finite meshes. Before any geometrically nonlinear implementation, verify residual directional derivatives, tangent symmetry for conservative loading, objectivity under large -rigid motion, zero stress-free geometric stiffness, and convergence of Newton -iterations. These checks are future-only and do not authorize a nonlinear procedure. +rigid motion, the first/second derivatives of `Phi`, the coordinate-map curvature +term, objective drilling, zero stress-free geometric stiffness, and convergence of +Newton iterations. These checks are future-only and do not authorize a nonlinear +procedure. ## 18. Numerical risks @@ -1404,17 +1690,18 @@ iterations. These checks are future-only and do not authorize a nonlinear proced | Membrane locking on distorted curved meshes | slow or nonuniform convergence | distortion/curvature sweeps; document original MITC4 limitation | | Volumetric locking | N/A for the approved plane-stress shell contract | do not infer a three-dimensional incompressible formulation | | Wrong tying pair or engineering-shear factor | swapped/incorrect shear and loss of patch consistency | component-level tying and patch tests | -| Reversed or degenerate Jacobian | invalid basis, sign, or energy | fail at every required point using scale-aware measures | +| Reversed or degenerate Jacobian | invalid basis, sign, or energy | common location inventory with `J/j_s/a_g/c_d/r_J/theta_w` | | Discontinuous shared director | artificial coupling or undefined frame | reject; require duplicate nodes at folds | | Fixed-axis tangent singularity | nondeterministic rotation transform | least-aligned global-axis construction | | Drilling coefficient too small | rank/conditioning failure | rank and conditioning sensitivity sweep | | Drilling coefficient too large | contaminated displacement/resultant | physical-output and separate-energy sensitivity sweep | -| Mixed-unit diagonal scale | unit-dependent stiffness | prohibit `min(all Kii)`; use rotationally compatible scale | +| Mixed-unit drilling or spectrum scale | unit-dependent stabilization/rank | use `rho_d,I` and `S^T K S`; prohibit raw mixed-unit comparison | | Misconstructed rigid test | false drill energy in a physical mode | use rigid translation plus tangent director change with `gamma=0` | | Recovery/stiffness mismatch | inconsistent energy and section output | same frames, tying, material, and thickness rule | | Location averaging | hidden sign/identity error | preserve exact location rows; no nodal extrapolation | | Treating S4R as reduced FESA integration | divergent kernel and unsupported hourglass behavior | one documented rule for S4/S4R | | Nonlinear director Hessian omission | inconsistent future tangent | exact first/second derivative of approved rotation map | +| Missing nonlinear 20-to-24 map Hessian | incomplete global tangent | require `Phi`, `A`, and coordinate-map curvature term | | Nonobjective nonlinear drill penalty | artificial large-rotation energy | future review of rotating frame and drill potential | ## 19. Evidence basis and applicability @@ -1461,41 +1748,65 @@ research brief remain the project source of truth. | Requirement group | Formulation coverage | Remaining owner | | --- | --- | --- | | `001-004`, `030`, `037` | linear-static boundary; S4/S4R one FESA path, source identity distinct | I/O, planning | -| `005`, `031-038` | 24 global DOFs; 20 physical plus four drill coordinates; transform, rank, energy | Numerical Review for drill choices | +| `005`, `031-038` | 24 global DOFs; 20 physical plus four drill coordinates; `rho_d,I`, scaled rank, energy | Numerical Review for drill calibration | | `006-010` | isotropic plane stress, one centered constant-thickness layer | I/O validation | -| `011-016` | connectivity-positive normals, deterministic area weighting, frames, Jacobian checks | Numerical Review for thresholds | -| `017-020` | global BC/CLOAD virtual work; parallel moment rejected; distributed loads excluded | I/O projection diagnostic | +| `011-016` | pairwise normals, deterministic averaging/frames, complete geometry measures | Numerical Review for calibrated thresholds | +| `017-020` | global 24-DOF virtual work; `rho_M <= 1e-12`; distributed loads excluded | I/O diagnostic/schema | | `024-029` | deterministic element buffers, partitioned linear lifecycle, full-residual reaction | planning | | `039-048` | nodal/global and shell/local recovery inventory, units, identities, energy split | I/O schema | -| `049-057` | invariants, patches, locking, geometry, curved shells, drilling, equilibrium | Numerical Review/reference/physics | +| `049-057` | scaled normalized invariants, patches, locking, geometry, curved shells, drilling, equilibrium | Numerical Review/reference/physics | | `058-064` | U blocking, UR warning-only comparison boundary | Numerical Review/reference model | | `065-072` | reference immutability and displacement-only evidence boundary acknowledged | reference model | +### 20.1 Numerical Review revision traceability + +The first review findings map to this revision as follows. + +`NR-C01` maps to Sections 9.2-9.3, 16.2, and 17.1/17.3. The geometry measures +are fixed while their thresholds still require calibration. + +`NR-C02` maps to Sections 12.2 and 17.4. The common `rho_d,I` coordinate and +candidate conversions are fixed while the stable plateau remains open. + +`NR-C03` maps to Sections 12.3/12.5 and 17.1. Element/global DOF scaling is fixed +while the numerical-rank threshold remains open. + +`NR-C04` maps to Sections 7.1-7.2, which fix the global 24-DOF weak form and +transpose notation. `NR-C05` maps to Sections 15.1-15.5 and 17.6, which expose the +physical-chart/global-map boundary and required map Hessian term. + +`NR-D01` maps to Sections 6.2 and 16.3, which fix the exact-zero case and +`rho_M <= 1e-12`. `NR-D02` maps to Sections 12.5 and 17.1, which fix the approved +normalized algebraic metrics. + ## 21. Open issues and downstream handoff ### 21.1 Blocking Numerical Review decisions 1. Select the dimensionally compatible drilling reference-scale family. -2. Select `alpha_d`, its sensitivity range, conditioning/rank acceptance, and the - physical-output contamination bound. +2. Select the nodewise `rho_d,I` stable plateau, nominal value, scaled + conditioning/rank acceptance, and physical-output contamination bound. 3. Define classification and warning behavior for `E_drill/E_phys`, including the zero/near-zero physical-energy case. 4. Approve `theta_smooth` after curved-mesh resolution sweeps; `20 degrees` is only the initial candidate. -5. Approve dimensionless Jacobian, distortion, and warpage measures and thresholds. -6. Approve the scale-aware tolerance used to reject a director-parallel nodal - moment. -7. Re-derive/check the exact `B_bar` component signs, tying-point inventory, - engineering-shear factors, rank, and bottom/top stress sign independently. +5. Calibrate thresholds for the fixed `J/j_s/a_g/c_d/r_J/theta_w` geometry + inventory using valid distortion/warp and collapsed negative sequences. +6. Jointly with Reference Model, approve the U mixed tolerance and nonblocking UR + large-error warning threshold. -Until these items are approved, the formulation is ready for Numerical Review but -not for Implementation Planning. +The first review already approved `rho_M <= 1e-12`, the normalized algebraic checks, +the MITC tying/component signs, constitutive law, quadrature, and recovery signs. +This revision is ready for Numerical Review rerun but not for Implementation +Planning until the six remaining evidence-backed decisions are closed. ### 21.2 I/O Definition handoff - Preserve source `S4`/`S4R` separately from internal `FESA-MITC4`. - Define exact keyword subset, section/material resolution, and fail-closed - diagnostics for director, folds, Jacobians, unsupported loads, and recovery. + diagnostics for director, folds, Jacobians, unsupported loads, and recovery; + encode exact-zero nodal moment separately and enforce `rho_M <= 1e-12` without + a denominator clamp. - Define exact HDF5 row schemas for global `U/UR`, `RF/RM`; four midsurface generalized-strain/resultant locations; bottom/middle/top stress positions; full residual/equilibrium; and separate energies. @@ -1522,8 +1833,9 @@ not for Implementation Planning. ### 21.5 Future nonlinear formulation handoff Before geometric-nonlinear implementation, separately approve the finite-director -chart and update law, objective drilling potential and scale update, exact -first/second director derivatives, nonlinear constitutive boundary, nodal-moment -work/load tangent, state ownership, convergence controls, and nonlinear output -contract. Section 15 supplies the residual/tangent structure but deliberately does -not close those product decisions. +chart and update law, the complete `Phi: R24 -> R20` map and its first/second +derivatives, objective drilling potential and scale update, exact first/second +director derivatives, nonlinear constitutive boundary, nodal-moment work/load +tangent, state ownership, convergence controls, and nonlinear output contract. +Section 15 supplies the physical-chart residual/tangent and the conditional global +mapping identity but deliberately does not close those product decisions.