baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

docs: revise MITC4 shell formulation after review

Browse Source
This commit is contained in:
KOKO\Mimi
2026-08-11 19:59:52 +09:00
parent 0a5aad427e
commit 6cb0302015
1 changed files with 439 additions and 127 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/formulations/mitc4-shell-formulation.md
+439 -127
View File
@@ -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_{e<f,\ e,f\in E(I)}\mathbf n_e\cdot\mathbf n_f>0.
$$
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.