diff --git a/docs/formulations/mitc4-shell-formulation.md b/docs/formulations/mitc4-shell-formulation.md new file mode 100644 index 0000000..3520521 --- /dev/null +++ b/docs/formulations/mitc4-shell-formulation.md @@ -0,0 +1,1529 @@ +# FESA Linear Static MITC4 Shell Formulation + +## Metadata + +- 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` +- status: `ready-for-numerical-review` +- owner_agent: `formulation-agent` +- date: `2026-08-11` +- product_execution_scope: `small-strain, small-rotation linear static only` +- future_formulation_scope: `geometrically nonlinear Total Lagrangian residual/tangent; not executable` + +## 1. Purpose, authority, and decision boundary + +This document is the mathematical and algorithmic contract for the FESA four-node +MITC4 shell element. It strengthens the earlier draft with the approved requirements +and research evidence. It does not define C++ APIs, parser classes, HDF5 paths, or a +reference-solver workflow. + +The normative product path is: + +- one linear-static step; +- small strain and small rotation; +- four-node bilinear degenerated-continuum shell geometry; +- 20 physical element coordinates, embedded in 24 global shell DOFs; +- MITC4 edge-midpoint assumed transverse shear; +- homogeneous isotropic linear elasticity, one centered layer, constant thickness; +- one documented `2 x 2 x 2` integration path; +- nodal `CLOAD` and prescribed global shell DOFs only. + +Section 15 records a future Total Lagrangian residual and consistent tangent. That +section is formulation evidence only. It does not expand the current analysis +procedure, mutable state, input subset, or implementation plan. + +Source labels `S4` and `S4R` both select this one FESA formulation by an approved +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. + +## 2. Scope and assumptions + +### 2.1 Included product formulation + +- Reference coordinates: `(xi, eta, zeta) in [-1,1]^3`. +- Source node order: `1=(-1,-1)`, `2=(+1,-1)`, `3=(+1,+1)`, + `4=(-1,+1)`. +- Global nodal DOF order: + `[UX, UY, UZ, URX, URY, URZ]`. +- Physical nodal fields: three midsurface translations and two director-tangent + rotations. +- Numerical nodal field: one director-parallel drilling coordinate. +- Material: finite `E > 0`, `-1 < nu < 0.5`. +- Section: one centered homogeneous isotropic layer with finite constant `t > 0`. +- Plane-stress condition `sigma_33 = 0`; no thickness stretch. +- User-consistent units. Angles are dimensionless. +- Initial directors are derived from source connectivity and smooth-patch averaging; + explicit normals are outside this feature. +- External loads are global nodal forces and non-drilling nodal moments. + +### 2.2 Excluded product formulation + +- Geometrically nonlinear execution, follower loads, dynamics, buckling, contact, + material nonlinearity, and thermal strain. +- Composite or multilayer sections, offsets, variable/distributed thickness, + orthotropy, and user material orientation. +- Pressure, gravity, body force, edge traction, `DLOAD`, or any distributed-load + input. +- MITC4+, MITC4/D, EAS/ANS variants other than the shear tying defined here, + reduced integration, hourglass control, and selective integration. +- A physical drilling strain, drilling stress/resultant, or drilling-direction + applied nodal moment. +- Shared-node sharp folds, hinges, shell intersections, and shell-beam drilling + transfer. A physical fold requires duplicate source nodes. + +### 2.3 Kinematic and constitutive assumptions + +The initial thickness director is a unit vector and thickness is a separate scalar. +Straight material lines through the thickness remain straight in the physical +five-DOF kinematics. The transverse normal strain and stress are excluded, while +transverse shear deformation is retained. The shear correction factor is + +$$ +kappa_s = \frac{5}{6}. +$$ + +The `5/6` value is restricted to this homogeneous rectangular-thickness contract. +It is not a general composite-shell rule. + +## 3. Coordinates, signs, and location identity + +### 3.1 Natural coordinates and positive face + +The positive thickness direction is `+zeta`. With the approved node order, the +positive element normal candidate follows + +$$ +\mathbf n_e \parallel \mathbf A_\xi \times \mathbf A_\eta. +$$ + +The bottom, middle, and top section positions are respectively + +$$ +z=-t/2,\qquad z=0,\qquad z=+t/2, +\qquad z=\frac{t}{2}\zeta. +$$ + +Positive section moment is defined by + +$$ +\mathbf M=\int_{-t/2}^{t/2} z\,\boldsymbol\sigma_m(z)\,dz. +$$ + +This convention fixes the curvature sign in Section 14. Reversing source +connectivity reverses the positive face and is not silently repaired. + +### 3.2 Integration-point local frame + +At each midsurface location, let the interpolated positive-thickness director be + +$$ +\overline{\mathbf d}=\frac{\sum_I N_I\mathbf d_I} +{\left\|\sum_I N_I\mathbf d_I\right\|}. +$$ + +Then define + +$$ +\mathbf e_3=\overline{\mathbf d}, +\qquad +\mathbf e_1= +\frac{\mathbf A_\xi-(\mathbf A_\xi\cdot\mathbf e_3)\mathbf e_3} + {\|\mathbf A_\xi-(\mathbf A_\xi\cdot\mathbf e_3)\mathbf e_3\|}, +$$ + +$$ +\mathbf e_2=\mathbf e_3\times\mathbf e_1. +$$ + +Thus local direction 1 follows the projected `+xi` covariant tangent, direction 2 +completes a right-handed frame, and direction 3 follows the interpolated positive +thickness director. The source-order surface normal must have positive dot product +with `e_3`; an opposed or degenerate pairing is invalid. +The frame is evaluated and stored separately at every required shell location; one +element-average frame is not substituted for warped geometry. Frame degeneracy is a +geometry error. + +All generalized strains, resultants, and section stresses in this document use this +local frame. Nodal `U/UR` and `RF/RM` remain global. + +### 3.3 Stable location identity + +No recovery values at different locations are averaged to hide mismatches. A shell +row is identified at least by source element identity, source element type, internal +formulation identity, in-plane natural coordinates, section position when +applicable, local frame/director, and component order. Exact HDF5 field names belong +to the I/O contract. + +## 4. Initial directors and nodal tangent frames + +### 4.1 Element normal candidate + +For element `e`, evaluate the midsurface covariant tangents at the element center: + +$$ +\mathbf A_\xi^e(0,0)=\sum_{I=1}^4 N_{I,\xi}(0,0)\mathbf X_I, +\qquad +\mathbf A_\eta^e(0,0)=\sum_{I=1}^4 N_{I,\eta}(0,0)\mathbf X_I. +$$ + +The source-order-positive candidate is + +$$ +\mathbf n_e= +\frac{\mathbf A_\xi^e\times\mathbf A_\eta^e} + {\|\mathbf A_\xi^e\times\mathbf A_\eta^e\|}. +$$ + +Its deterministic surface-area weight is + +$$ +A_e=\sum_{g=1}^{4}w_g +\left\|\mathbf A_\xi(\xi_g,\eta_g) +\times\mathbf A_\eta(\xi_g,\eta_g)\right\|, +$$ + +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 + +$$ +\mathbf s_I=\sum_{e\in E(I)} A_e\mathbf n_e, +\qquad +\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 + +$$ +\theta_{eI}=\cos^{-1}\!\left( +\operatorname{clamp}(\mathbf n_e\cdot\mathbf d_I,-1,1)\right) +$$ + +must satisfy an approved smooth-patch bound `theta_smooth`. The research value +`20 degrees` is only the first Numerical Review candidate. It is not fixed here. +Zero or near-zero cross products and averaged vectors fail using approved +scale-aware tolerances. They are never replaced with a global axis. + +### 4.3 Deterministic nodal tangent frame + +Let `(g_x,g_y,g_z)` be the global Cartesian basis. Choose the least-aligned basis +axis + +$$ +\mathbf g_*= +\underset{\mathbf g\in\{\mathbf g_x,\mathbf g_y,\mathbf g_z\}} +{\arg\min}\ |\mathbf g\cdot\mathbf d_I|, +$$ + +with ties resolved in the fixed order `x`, `y`, `z`. Then + +$$ +\mathbf a_I= +\frac{\mathbf g_*-(\mathbf g_*\cdot\mathbf d_I)\mathbf d_I} + {\|\mathbf g_*-(\mathbf g_*\cdot\mathbf d_I)\mathbf d_I\|}, +\qquad +\mathbf b_I=\mathbf d_I\times\mathbf a_I. +$$ + +The nodal frame `(a_I,b_I,d_I)` is orthonormal and right-handed: +`a_I x b_I=d_I`. + +## 5. Primary variables and six-to-five DOF embedding + +### 5.1 Global and local nodal coordinates + +At node `I`, the global shell vector is + +$$ +\mathbf q_I^g= +\begin{bmatrix} +\mathbf u_I\\ \boldsymbol\theta_I^g +\end{bmatrix} += +\begin{bmatrix} +U_X&U_Y&U_Z&UR_X&UR_Y&UR_Z +\end{bmatrix}^T. +$$ + +Define + +$$ +\mathbf R_I= +\begin{bmatrix}\mathbf a_I&\mathbf b_I&\mathbf d_I\end{bmatrix}, +\qquad +\begin{bmatrix}\alpha_I\\\beta_I\\\gamma_I\end{bmatrix} +=\mathbf R_I^T\boldsymbol\theta_I^g. +$$ + +`alpha` and `beta` are physical director-tangent rotations. `gamma` is the numerical +drilling coordinate. Under the small-rotation convention + +$$ +\delta\mathbf d_I +=\boldsymbol\theta_I^g\times\mathbf d_I +=\beta_I\mathbf a_I-\alpha_I\mathbf b_I, +$$ + +and the `gamma_I d_I` part vanishes identically from director change. + +### 5.2 Element transforms + +Collect the four global nodal vectors as `q_g in R^24`. Define the block transform + +$$ +\mathbf T=\operatorname{blockdiag}(\mathbf T_1,\ldots,\mathbf T_4), +\qquad +\mathbf T_I= +\begin{bmatrix} +\mathbf I_3&\mathbf 0\\ +\mathbf 0&\mathbf R_I^T +\end{bmatrix}. +$$ + +Let `P in R^(20 x 24)` select, for each node, +`[u_X,u_Y,u_Z,alpha,beta]`, and let `D_gamma in R^(4 x 24)` select the four +`gamma` entries. Then + +$$ +\boxed{\mathbf q_{20}=\mathbf T_p\mathbf q_g, +\qquad \mathbf T_p=\mathbf P\mathbf T}, +$$ + +$$ +\boxed{\boldsymbol\gamma=\mathbf T_d\mathbf q_g, +\qquad \mathbf T_d=\mathbf D_\gamma\mathbf T}. +$$ + +`T_p` is `20 x 24` and `T_d` is `4 x 24`. The physical strain operator acts only +on `q_20`; `gamma` must not enter physical strain or recovery. + +### 5.3 Virtual work preservation + +For any physical local force `f_20` and drill force `m_d`, + +$$ +\delta\mathbf q_{20}^T\mathbf f_{20} ++\delta\boldsymbol\gamma^T\mathbf m_d +=\delta\mathbf q_g^T +\left(\mathbf T_p^T\mathbf f_{20}+\mathbf T_d^T\mathbf m_d\right). +$$ + +This identity is the required transformation-energy check. A physical rigid mode is +constructed from rigid midsurface translations and tangent-plane director changes; +blindly assigning the complete spatial rigid-rotation vector, including its normal +projection, to every `UR` entry would excite the numerical drill coordinate and is +not the five-DOF rigid mode. + +## 6. Strong form and boundary conditions + +### 6.1 Three-dimensional parent statement + +On the reference shell volume `Omega_0`, the linear elastic parent problem is + +$$ +\nabla\cdot\boldsymbol\sigma=\mathbf0 +\quad\text{in }\Omega_0, +\qquad +\boldsymbol\sigma=\mathbb C_{ps}:\boldsymbol\epsilon, +$$ + +with + +$$ +\boldsymbol\epsilon=\operatorname{sym}(\nabla\mathbf u), +\qquad \sigma_{33}=0. +$$ + +The continuum strong form explains the element lineage. In the approved product +path there is no body force and no surface/edge traction input. + +### 6.2 Essential and natural data + +Prescribed global shell DOFs form the Dirichlet data. The only external generalized +forces are aggregated nodal `CLOAD` values. For a nodal moment `M_I`, define + +$$ +M_{d,I}=\mathbf d_I\cdot\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. + +Distributed natural-boundary terms may appear in the parent continuum derivation, +but their existence is not an executable load contract. + +## 7. Weak form + +### 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 + +$$ +\delta W_{int}^{phys}-\delta W_{ext}=0 +\qquad\forall\delta\mathbf q_{20}\in V_{20}^0, +$$ + +where + +$$ +\delta W_{int}^{phys} +=\int_{\Omega_0} +\delta\overline{\boldsymbol\epsilon}^{,T} +\mathbf C_5 +\overline{\boldsymbol\epsilon}\,dV, +$$ + +and `overbar` denotes the MITC-projected shear field. In the current product path + +$$ +\delta W_{ext}=\delta\mathbf q_g^T\mathbf f_e^{CLOAD}. +$$ + +### 7.2 Numerical drilling potential + +Drilling regularization is algebraically separate: + +$$ +\Pi_d=\frac12\boldsymbol\gamma^T\mathbf K_d^l\boldsymbol\gamma, +\qquad +\delta\Pi_d=\delta\boldsymbol\gamma^T +\mathbf K_d^l\boldsymbol\gamma. +$$ + +It does not approximate a shell strain energy and produces no physical generalized +strain, resultant, or stress. + +## 8. Bilinear discretization and degenerated geometry + +### 8.1 Shape functions + +With nodal signs `(xi_I,eta_I)`, + +$$ +N_I(\xi,\eta)=\frac14(1+\xi_I\xi)(1+\eta_I\eta), +$$ + +or explicitly + +$$ +\begin{aligned} +N_1&=\tfrac14(1-\xi)(1-\eta),& +N_2&=\tfrac14(1+\xi)(1-\eta),\\ +N_3&=\tfrac14(1+\xi)(1+\eta),& +N_4&=\tfrac14(1-\xi)(1+\eta). +\end{aligned} +$$ + +They must satisfy + +$$ +\sum_I N_I=1, +\qquad N_I(\xi_J,\eta_J)=\delta_{IJ}, +\qquad \sum_I N_{I,\xi}=\sum_I N_{I,\eta}=0. +$$ + +### 8.2 Reference geometry + +For element thickness `t`, + +$$ +\boxed{ +\mathbf X(\xi,\eta,\zeta) +=\sum_{I=1}^4N_I\mathbf X_I ++\frac{t\zeta}{2}\sum_{I=1}^4N_I\mathbf d_I +}. +$$ + +The midsurface is `X_0(xi,eta)=X(xi,eta,0)`. The initial covariant bases are + +$$ +\mathbf G_\alpha=\frac{\partial\mathbf X}{\partial\xi^\alpha}, +\qquad (\xi^1,\xi^2,\xi^3)=(\xi,\eta,\zeta). +$$ + +### 8.3 Linear displacement interpolation + +Using `delta d_I=beta_I a_I-alpha_I b_I`, + +$$ +\boxed{ +\mathbf u(\xi,\eta,\zeta) +=\sum_{I=1}^4N_I\mathbf u_I ++\frac{t\zeta}{2}\sum_{I=1}^4N_I +(\beta_I\mathbf a_I-\alpha_I\mathbf b_I) +}. +$$ + +Therefore + +$$ +\mathbf u=\mathbf H(\xi,\eta,\zeta)\mathbf q_{20}, +\qquad \mathbf H\in\mathbb R^{3\times20}. +$$ + +This interpolation reproduces constant translation. Its physical rigid-rotation +test must combine `u_I=omega x X_I` with +`delta d_I=omega x d_I` projected into `(a_I,b_I)`. The corresponding global nodal +rotation variable is the tangent projection + +$$ +\boldsymbol\theta_I^g +=\boldsymbol\omega-(\boldsymbol\omega\cdot\mathbf d_I)\mathbf d_I, +$$ + +so that `gamma_I=0`. + +## 9. Mapping, Jacobian, and geometry validity + +### 9.1 Covariant and contravariant bases + +Define + +$$ +\mathbf J= +\begin{bmatrix}\mathbf G_\xi&\mathbf G_\eta&\mathbf G_\zeta\end{bmatrix}, +\qquad J=\det\mathbf J, +$$ + +and contravariant bases by + +$$ +\mathbf G^\alpha\cdot\mathbf G_\beta=\delta^\alpha_\beta. +$$ + +The volume measure is + +$$ +dV=J\,d\xi\,d\eta\,d\zeta. +$$ + +### 9.2 Required validation locations + +`J` and all basis quantities must be finite and valid at: + +- all eight `2 x 2 x 2` stiffness points; +- the four midsurface MITC tying locations; +- any additional bottom/middle/top recovery evaluation location before recovery is + committed. + +Duplicate nodes, self-intersection, degenerate midsurface area, and reversed +connectivity are separate fail-closed geometry errors. + +### 9.3 Scale-aware measures + +A dimensionless local determinant measure is + +$$ +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. + +## 10. Linear kinematics and MITC4 shear projection + +### 10.1 Direct covariant strain + +For a physical element coordinate `q_A`, let + +$$ +\mathbf h_A=\frac{\partial\mathbf u}{\partial q_A}, +\qquad +\mathbf h_{A,\alpha}=\frac{\partial\mathbf h_A}{\partial\xi^\alpha}. +$$ + +The direct covariant small-strain column is + +$$ +B^{DI}_{\alpha\beta,A} +=\frac12\left( +\mathbf G_\alpha\cdot\mathbf h_{A,\beta} ++\mathbf G_\beta\cdot\mathbf h_{A,\alpha} +\right). +$$ + +Equivalently, + +$$ +\boldsymbol\epsilon^{DI} +=\operatorname{sym}\left( +\sum_{\alpha=1}^3 +\mathbf u_{,\alpha}\otimes\mathbf G^\alpha +\right). +$$ + +The transverse normal component is excluded by the shell constitutive contract. + +### 10.2 Canonical tying locations + +Use + +$$ +T_{\xi-}=(0,-1,0),\qquad T_{\xi+}=(0,+1,0), +$$ + +$$ +T_{\eta-}=(-1,0,0),\qquad T_{\eta+}=(+1,0,0). +$$ + +The first pair supplies the covariant `xi-zeta` shear and the second pair supplies +the covariant `eta-zeta` shear. Labels such as A/B/C/D are avoided because their +assignment varies with published node conventions. + +### 10.3 Assumed transverse shear + +At any in-plane point, + +$$ +\boxed{ +\overline\epsilon_{\xi\zeta}(\xi,\eta) +=\frac{1-\eta}{2}\epsilon^{DI}_{\xi\zeta}(T_{\xi-}) ++\frac{1+\eta}{2}\epsilon^{DI}_{\xi\zeta}(T_{\xi+}) +}, +$$ + +$$ +\boxed{ +\overline\epsilon_{\eta\zeta}(\xi,\eta) +=\frac{1-\xi}{2}\epsilon^{DI}_{\eta\zeta}(T_{\eta-}) ++\frac{1+\xi}{2}\epsilon^{DI}_{\eta\zeta}(T_{\eta+}) +}. +$$ + +All in-plane covariant components remain direct. The assumed covariant tensor is +reconstructed with the contravariant bases: + +$$ +\overline{\boldsymbol\epsilon} +=\sum_{\alpha,\beta} +\overline\epsilon_{\alpha\beta} +\mathbf G^\alpha\otimes\mathbf G^\beta. +$$ + +### 10.4 Local engineering strain vector + +Transform the tensor to `(e_1,e_2,e_3)` and define + +$$ +\boxed{ +\overline{\mathbf e}= +\begin{bmatrix} +\epsilon_{11}&\epsilon_{22}&\gamma_{12}&\gamma_{13}&\gamma_{23} +\end{bmatrix}^T +=\overline{\mathbf B}\mathbf q_{20}}, +$$ + +where `gamma_ij=2 epsilon_ij` and +`B_bar in R^(5 x 20)`. The normative construction of `B_bar` is: + +1. form each direct covariant strain column; +2. replace only its two transverse-shear covariant components by the same tying + projection; +3. reconstruct the Cartesian tensor with the contravariant bases; +4. project it into the stored local frame and apply engineering-shear factors. + +Strain, first variation, and future second variation must use the same linear MITC +projection. Using tied shear in residual but direct shear in stiffness is forbidden. + +## 11. Constitutive contract + +### 11.1 Plane-stress material matrix + +Let + +$$ +G=\frac{E}{2(1+\nu)}, +\qquad +\mathbf C_{ps}=\frac{E}{1-\nu^2} +\begin{bmatrix} +1&\nu&0\\ +\nu&1&0\\ +0&0&(1-\nu)/2 +\end{bmatrix}. +$$ + +For the engineering order of Section 10, + +$$ +\boxed{ +\mathbf C_5= +\begin{bmatrix} +\mathbf C_{ps}&\mathbf0\\ +\mathbf0&\kappa_sG\mathbf I_2 +\end{bmatrix}}, +$$ + +$$ +\mathbf s= +\begin{bmatrix} +\sigma_{11}&\sigma_{22}&\tau_{12}&\tau_{13}&\tau_{23} +\end{bmatrix}^T +=\mathbf C_5\overline{\mathbf e}. +$$ + +`C_5` is symmetric positive definite for the approved material range. `sigma_33=0` +is an assumption, not an emitted stress row. + +There are no material state variables, history variables, or constitutive update in +the approved linear elastic path. + +### 11.2 Centered-section matrices as a cross-check + +For an exactly linear in-plane strain through thickness, + +$$ +\mathbf A=t\mathbf C_{ps},\qquad +\mathbf B=\mathbf0,\qquad +\mathbf D=\frac{t^3}{12}\mathbf C_{ps},\qquad +\mathbf A_s=\kappa_sGt\mathbf I_2. +$$ + +These matrices are recovery and verification identities. The normative element +stiffness remains the source-faithful volume integration in Section 12, so curved or +warped mapping terms are not silently replaced by a flat analytical section kernel. + +## 12. Linear element equations and drilling stabilization + +### 12.1 Physical 20-DOF kernel + +Using the quadrature of Section 13, + +$$ +\boxed{ +\mathbf K_{20} +=\int_{\Omega_0}\overline{\mathbf B}^T +\mathbf C_5\overline{\mathbf B}\,dV +\approx\sum_{g=1}^{8} +\overline{\mathbf B}_g^T\mathbf C_5\overline{\mathbf B}_g +J_gw_g}, +$$ + +$$ +\mathbf K_{20}\in\mathbb R^{20\times20}, +\qquad +\mathbf f_{int}^{20}=\mathbf K_{20}\mathbf q_{20}. +$$ + +The physical global contribution is + +$$ +\boxed{ +\mathbf K_{phys}^{24}=\mathbf T_p^T\mathbf K_{20}\mathbf T_p}, +\qquad +\mathbf f_{phys}^{24}=\mathbf K_{phys}^{24}\mathbf q_g. +$$ + +### 12.2 Drilling candidate contract + +Let + +$$ +\mathbf K_d^l=\operatorname{diag}(k_{d,1},k_{d,2},k_{d,3},k_{d,4}), +\qquad k_{d,I}>0, +$$ + +with rotational-stiffness units `force*length`. Write + +$$ +k_{d,I}=\alpha_d k_{ref,I}, +$$ + +where both the dimensionless coefficient `alpha_d` and the reference family remain +unapproved until Numerical Review. Dimensionally compatible candidates carried from +research are: + +1. transverse-shear/area transition family + +$$ +k_{ref,I}^{(A)}= +\frac{GtA_{eI}}{1+qA_{eI}/t^2}, +\qquad q=2.5\times10^{-5}, +$$ + +where + +$$ +A_{eI}=\int_{A_e}N_I\,dA +\approx\sum_{g=1}^{4}N_I(\xi_g,\eta_g) +\|\mathbf A_\xi\times\mathbf A_\eta\|_g w_g; +$$ + +2. isotropic bending rigidity + +$$ +k_{ref}^{(B)}=D_{iso}=\frac{Et^3}{12(1-\nu^2)}; +$$ + +3. a documented positive statistic formed only from the physical rotational block + of `K_20`, whose entries all have `force*length` units. + +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. + +### 12.3 Stabilized 24-DOF matrix + +The global drilling contribution and total element matrix are + +$$ +\boxed{ +\mathbf K_{drill}^{24}=\mathbf T_d^T\mathbf K_d^l\mathbf T_d}, +$$ + +$$ +\boxed{ +\mathbf K_e^{24}=\mathbf K_{phys}^{24}+\mathbf K_{drill}^{24}}, +$$ + +$$ +\mathbf f_{int}^e=\mathbf K_e^{24}\mathbf q_g, +\qquad +\mathbf r_e=\mathbf f_{int}^e-\mathbf f_e^{CLOAD}. +$$ + +No element distributed-load integral exists in the approved product routine. +Mass and damping matrices are `N/A` for this linear-static feature. + +For a valid free isolated element, the expected physical rank is 14. Embedding it in +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. + +### 12.4 Energy split + +The element energies are + +$$ +\boxed{ +E_{phys}^e=\frac12\mathbf q_g^T +\mathbf K_{phys}^{24}\mathbf q_g +=\frac12\mathbf q_{20}^T\mathbf K_{20}\mathbf q_{20}}, +$$ + +$$ +\boxed{ +E_{drill}^e=\frac12\mathbf q_g^T +\mathbf K_{drill}^{24}\mathbf q_g +=\frac12\boldsymbol\gamma^T\mathbf K_d^l\boldsymbol\gamma}. +$$ + +Both have units `force*length` and are aggregated separately in stable source order. +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. + +## 13. Numerical integration + +### 13.1 Stiffness rule + +Use tensor-product Gauss points + +$$ +\xi,\eta,\zeta\in\left\{-\frac1{\sqrt3},+\frac1{\sqrt3}\right\}, +\qquad w_\xi=w_\eta=w_\zeta=1. +$$ + +Thus the physical kernel has four midsurface locations and two thickness locations, +for eight volume evaluations. At each in-plane Gauss location: + +1. evaluate and validate the midsurface local frame; +2. evaluate the four tying values once at `zeta=0`; +3. use the same tied covariant shear interpolation at both thickness points; +4. evaluate the remaining strain, mapping, material, and volume measure at the + actual thickness point. + +`S4` and `S4R` use this identical rule. One-point reduced integration and Abaqus +hourglass control are not selected by the `S4R` source label. + +### 13.2 Prohibited quadrature substitutions + +- replacing `2 x 2` midsurface integration with `1 x 1`; +- using direct shear at any stiffness point; +- evaluating tying values from a different geometry or director inventory; +- averaging failed or missing tying values; +- changing the rule based on source `S4` versus `S4R`; +- using higher-order quadrature in production without reopening Formulation and + Numerical Review. + +Higher-order or analytical integration may be used only as an independent +verification cross-check. + +## 14. Output recovery + +### 14.1 Nodal solution, reaction, and equilibrium + +Store at every source node, in global coordinates, + +$$ +[U1,U2,U3,UR1,UR2,UR3] +$$ + +and + +$$ +[RF1,RF2,RF3,RM1,RM2,RM3]. +$$ + +The reaction/equilibrium vector is the assembled full residual + +$$ +\mathbf r=\mathbf K\mathbf d-\mathbf F. +$$ + +Constrained components are physical reactions; free components are residual +evidence. Reaction is not reconstructed by summing section resultants. + +For a fixed global reference point `X_o`, the force and moment balance evidence is + +$$ +\mathbf e_F=\sum_I\left(\mathbf F_I^{CLOAD}+\mathbf R_I\right), +$$ + +$$ +\mathbf e_M=\sum_I\left[ +(\mathbf X_I-\mathbf X_o)\times +(\mathbf F_I^{CLOAD}+\mathbf R_I) ++\mathbf M_I^{CLOAD}+\mathbf{RM}_I\right]. +$$ + +Only constrained residual components contribute to `R/RM`; free components remain +separate residual evidence. Normalization and acceptance thresholds belong to +Numerical Review and the I/O verification contract. + +### 14.2 Generalized strain by thickness moments + +At each of the four midsurface Gauss locations, let + +$$ +\mathbf e_m(z)= +\begin{bmatrix}\epsilon_{11}&\epsilon_{22}&\gamma_{12}\end{bmatrix}^T, +\qquad +\boldsymbol\gamma_s(z)= +\begin{bmatrix}\gamma_{13}&\gamma_{23}\end{bmatrix}^T, +$$ + +evaluated in the same stored local frame. Define + +$$ +\boldsymbol\epsilon_0= +\frac1t\int_{-t/2}^{t/2}\mathbf e_m(z)\,dz, +$$ + +$$ +\boldsymbol\kappa= +\frac{12}{t^3}\int_{-t/2}^{t/2}z\,\mathbf e_m(z)\,dz, +$$ + +$$ +\boldsymbol\gamma_0= +\frac1t\int_{-t/2}^{t/2}\boldsymbol\gamma_s(z)\,dz. +$$ + +Use the same two-point thickness quadrature: + +$$ +\int_{-t/2}^{t/2}f(z)\,dz +\approx\frac t2\sum_{h=1}^2w_h f(t\zeta_h/2). +$$ + +The mandatory generalized strain order is + +$$ +\boxed{ +[E11,E22,G12,K11,K22,K12,G13,G23] +=[\boldsymbol\epsilon_0,\boldsymbol\kappa, +\boldsymbol\gamma_0]}. +$$ + +Membrane and shear strains are dimensionless; curvature has units `1/length`. +For an exactly linear in-plane field, +`e_m(z)=epsilon_0+z kappa`, fixing the sign convention. + +### 14.3 Section resultants + +With + +$$ +\boldsymbol\sigma_m= +\begin{bmatrix}\sigma_{11}&\sigma_{22}&\tau_{12}\end{bmatrix}^T, +\qquad +\boldsymbol\tau_s= +\begin{bmatrix}\tau_{13}&\tau_{23}\end{bmatrix}^T, +$$ + +define + +$$ +\mathbf N=\int_{-t/2}^{t/2}\boldsymbol\sigma_m\,dz, +\qquad +\mathbf M=\int_{-t/2}^{t/2}z\boldsymbol\sigma_m\,dz, +\qquad +\mathbf Q=\int_{-t/2}^{t/2}\boldsymbol\tau_s\,dz. +$$ + +The mandatory order is + +$$ +\boxed{ +[N11,N22,N12,M11,M22,M12,Q13,Q23] +=[\mathbf N,\mathbf M,\mathbf Q]}. +$$ + +`N` and `Q` have units `force/length`; `M` has units `force`. For the centered +homogeneous linear-through-thickness case, + +$$ +\mathbf N=\mathbf A\boldsymbol\epsilon_0, +\qquad +\mathbf M=\mathbf D\boldsymbol\kappa, +\qquad +\mathbf Q=\mathbf A_s\boldsymbol\gamma_0, +$$ + +which is a required recovery cross-check. + +### 14.4 Bottom, middle, and top in-plane stress + +At each midsurface Gauss location, perform additional recovery evaluations at +`zeta=-1,0,+1` and emit only + +$$ +\boxed{[S11,S22,S12]}. +$$ + +These are direct section-position evaluations, not extrapolations from a different +location. `S33=0` is documented but not emitted. `S13/S23` point stresses are not +emitted because the approved external contract carries transverse shear through +`Q13/Q23`. + +### 14.5 No nodal extrapolation + +Element generalized quantities and section stresses remain at their defined shell +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. + +### 15.1 Total Lagrangian kinematics + +Let the physical nonlinear coordinate vector retain three translations and two +director parameters per node. In the reference configuration, + +$$ +\mathbf x(\xi,\eta,\zeta;\mathbf q) +=\sum_I N_I(\mathbf X_I+\mathbf u_I) ++\frac{t\zeta}{2}\sum_I N_I\mathbf d_I(\boldsymbol\phi_I), +$$ + +where `d_I(phi_I)` is a unit-director update defined by an approved finite-rotation +map. A candidate two-parameter exponential chart is + +$$ +\mathbf d_I(\boldsymbol\phi_I)= +\exp\!\left([ +\phi_{1I}\mathbf a_I+\phi_{2I}\mathbf b_I]_{\times}\right) +\mathbf d_I^0. +$$ + +The deformation gradient and Green-Lagrange strain are + +$$ +\mathbf F=\mathbf x_{,\alpha}\otimes\mathbf G^\alpha, +\qquad +\mathbf E=\frac12(\mathbf F^T\mathbf F-\mathbf I). +$$ + +For physical coordinates `q_A` and `q_B`, define + +$$ +\mathbf h_A=\frac{\partial\mathbf x}{\partial q_A}, +\qquad +\mathbf h_{AB}=\frac{\partial^2\mathbf x} +{\partial q_A\partial q_B}. +$$ + +The covariant Green-Lagrange components and their exact first and second derivatives +are + +$$ +E_{\alpha\beta}=\frac12\left( +\mathbf x_{,\alpha}\cdot\mathbf x_{,\beta} +-\mathbf G_\alpha\cdot\mathbf G_\beta\right), +$$ + +$$ +\boxed{ +E_{\alpha\beta,A}=\frac12\left( +\mathbf h_{A,\alpha}\cdot\mathbf x_{,\beta} ++\mathbf x_{,\alpha}\cdot\mathbf h_{A,\beta}\right)}, +$$ + +$$ +\boxed{ +\begin{aligned} +E_{\alpha\beta,AB}=\frac12\big(& +\mathbf h_{AB,\alpha}\cdot\mathbf x_{,\beta} ++\mathbf x_{,\alpha}\cdot\mathbf h_{AB,\beta}\\ +&+\mathbf h_{A,\alpha}\cdot\mathbf h_{B,\beta} ++\mathbf h_{B,\alpha}\cdot\mathbf h_{A,\beta} +\big). +\end{aligned}} +$$ + +For translational coordinates, `h_AB=0`. Rotation-rotation blocks retain the exact +second derivative of the director map; mixed-node blocks of `h_AB` are zero for the +nodal interpolation above. The quadratic products of `h_A` and `h_B` remain and form +the stress-dependent geometric contribution even when `h_AB=0`. + +Apply the same covariant MITC projection to the two transverse shear components of +`E`, its first variation, and its second variation: + +$$ +\overline{\mathbf E}=\mathcal P_{MITC}(\mathbf E), +$$ + +$$ +\delta\overline{\mathbf E} +=\mathcal P_{MITC}(\delta\mathbf E), +\qquad +\delta\Delta\overline{\mathbf E} +=\mathcal P_{MITC}(\delta\Delta\mathbf E). +$$ + +The exact first and second derivatives of the chosen director map are mandatory. +The mixed director Hessian is not assumed zero merely because the linearized +small-rotation expression is linear. + +Tensor-to-engineering transformation for the future material law uses the fixed +initial local frame at each reference integration point. Consequently `B_A` and +`G_AB` below are obtained by applying, in order, the MITC projection, the fixed +reference-frame tensor transformation, and the engineering-shear factors to the +explicit derivatives above. A future rotating output frame is a recovery decision +and must not be inserted into the constitutive tangent without its own derivatives. + +### 15.2 Residual + +Use second Piola-Kirchhoff stress `S` work-conjugate to `E`. In engineering-vector +notation, define + +$$ +\mathbf B_A(\mathbf q)= +\frac{\partial\overline{\mathbf e}}{\partial q_A}. +$$ + +For the geometrically nonlinear but materially linear candidate, + +$$ +\mathbf s=\mathbf C_5\overline{\mathbf e}, +\qquad +\mathbf C_T=\mathbf C_5, +$$ + +interpreted as a second-Piola/Green-Lagrange St. Venant-Kirchhoff-type shell law in +the fixed reference frame with the same plane-stress and shear-correction boundary +as the linear kernel. A different finite-strain or material-nonlinear law requires a +new constitutive formulation decision. + +Then the physical internal residual is + +$$ +\boxed{ +r_A^{phys}(\mathbf q) +=\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). +$$ + +### 15.3 Consistent material and geometric tangent + +Let + +$$ +\mathbf C_T=\frac{\partial\mathbf s} +{\partial\overline{\mathbf e}}, +\qquad +\mathbf G_{AB}=\frac{\partial^2\overline{\mathbf e}} +{\partial q_A\partial q_B}. +$$ + +The consistent physical tangent is + +$$ +\boxed{ +K_{AB}^{phys}=K_{AB}^{mat}+K_{AB}^{geo}}, +$$ + +$$ +\boxed{ +K_{AB}^{mat}=\int_{\Omega_0} +\mathbf B_A^T\mathbf C_T\mathbf B_B\,dV}, +$$ + +$$ +\boxed{ +K_{AB}^{geo}=\int_{\Omega_0} +\mathbf s^T\mathbf G_{AB}\,dV}. +$$ + +`K_geo` contains the stress-dependent second variation of Green-Lagrange strain, +including consistent finite-director derivatives. The future total tangent is + +$$ +\boxed{ +\mathbf K_T= +\mathbf K^{mat}+\mathbf K^{geo}+\mathbf K^{drill}-\mathbf K^{load}}. +$$ + +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, + +$$ +\Pi_d=\frac12\boldsymbol\gamma(\mathbf q)^T +\mathbf K_d^l(\mathbf q)\boldsymbol\gamma(\mathbf q), +\qquad +\mathbf r^{drill}=\frac{\partial\Pi_d}{\partial\mathbf q}, +\qquad +\mathbf K^{drill}=\frac{\partial^2\Pi_d}{\partial\mathbf q^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. + +### 15.5 Newton equation and linear limit + +At iteration `k`, + +$$ +\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. +$$ + +At the undeformed, stress-free state with the small-rotation director derivative, + +$$ +\mathbf K^{geo}=\mathbf0, +$$ + +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`. + +## 16. Algorithm pseudocode + +### 16.1 Deterministic preprocessing + +```text +for each shell element in stable source order: + validate four distinct source nodes and source-order geometry + 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 + 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) +``` + +### 16.2 Linear element kernel + +```text +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 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 + 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] + K20 += B_bar^T * C5 * B_bar * J * weight + +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 +Kdrill24 = T_d^T * Kd_local * T_d +Ke24 = Kphys24 + Kdrill24 +fint24 = Ke24 * q_g +residual24 = fint24 - f_CLOAD +return matrices, residual, transforms, frames, and separate energy operators +``` + +### 16.3 Global linear-static lifecycle + +```text +assemble all Ke24 contributions with stable element-local COO ordering +partition full K into Kff, Kfc, Kcf, Kcc in stable free/constrained order +factorize Kff before load assembly +assemble and deterministically aggregate nodal CLOAD +solve Kff * df = Ff - Kfc * dc +reconstruct full displacement d +compute full residual r = K*d - F +recover shell rows and physical/drilling energies in stable source order +validate complete finite candidate state/output, then commit +``` + +### 16.4 Recovery + +```text +for each element and each 2x2 midsurface location in fixed order: + evaluate two thickness Gauss strains/stresses with the stiffness formulation + integrate thickness zeroth/first moments for generalized strain/resultants + independently evaluate S11,S22,S12 at zeta = -1, 0, +1 + attach exact natural coordinates, section position, frame, and source identity + +recover nodal global U/UR and full-residual RF/RM +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) + for several normalized perturbation directions p and decreasing h: + compare K_T*p with [r(q+h*p)-r(q-h*p)]/(2*h) + require the error to decrease in the expected truncation range +``` + +## 17. Verification contract + +### 17.1 Algebraic and geometry checks + +- 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. +- Transformation work/energy invariance. +- Physical rigid modes satisfy normalized `Kphys24 r=0`. +- 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. + +### 17.2 Patch and sign checks + +Independently verify: + +- constant membrane strain/stress and `N` sign; +- pure bending about both local axes, `K11/K22`, `M`, and bottom/top stress sign; +- constant transverse shear and `Q13/Q23` order; +- pure twist and `K12/M12` convention; +- zero strain/resultant/stress contribution from a pure drilling vector. + +### 17.3 Locking, distortion, and curved shells + +- Thin and thick plate/shell mesh and thickness sequences are required; one + displacement on one mesh is insufficient. +- Distorted and warped valid quadrilaterals must be swept through approved geometry + measures. +- Original MITC4 controls transverse-shear locking but can retain membrane locking + in distorted curved meshes. This is a known limitation, not permission to add + MITC4+. +- Preferred nodal-load-compatible curved benchmarks are the pinched cylinder and + NAFEMS LE3 hemispherical shell. Scordelis-Lo is admissible only after an equivalent + nodal-load adaptation is documented. + +### 17.4 Drilling sensitivity + +For every candidate reference scale, sweep dimensionless coefficients around the +research order of magnitude and record: + +- free-element rank and 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. + +### 17.5 Reference-comparison boundary + +Abaqus comparisons block only on matched global `U1/U2/U3` rows under the approved +mixed tolerance decided downstream. `UR1/UR2/UR3` is fully reported and may emit a +deterministic nonblocking large-error warning, but it does not change pass/fail. +FESA `S4` and `S4R` inputs must produce the same internal numerical rows for identical +supported models while preserving distinct source metadata. Abaqus S4 and S4R are +not expected to be numerically identical on finite meshes. + +### 17.6 Future nonlinear tangent verification + +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. + +## 18. Numerical risks + +| Risk | Consequence | Required control | +| --- | --- | --- | +| Transverse-shear locking | overly stiff thin-shell response | exact edge-midpoint MITC projection and thickness/mesh convergence | +| 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 | +| 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 | +| 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 | +| Nonobjective nonlinear drill penalty | artificial large-rotation energy | future review of rotating frame and drill potential | + +## 19. Evidence basis and applicability + +### 19.1 Repository research sources + +The detailed source tiers, extracted facts, page references, benchmark provenance, +and evidence limits are owned by +`docs/research/linear-static-mitc4-shell-research.md`. The primary local source set +under `docs/reference-papers/MITC4/` includes: + +- `AContinuumMechanicsBasedFourNodeShell_001.md` and `_002.md`; +- `FourNodeQuadrilateralShellElementMITC4_001.md`; +- `MITC공부_001.md` and `_002.md`; +- `유한요소해석법을이용한쉘구조물의동적좌굴해석_001.md` through `_008.md`; +- `쉘구조물의유한요소해석에대하여_001.md` and `_002.md`. + +The evidence-backed core is the degenerated-continuum geometry, 20 physical DOFs, +edge-midpoint covariant shear tying, plane-stress constitutive response, and +`2 x 2` midsurface plus two-point thickness integration. Six-global-DOF drilling +regularization is a FESA interface/stability layer and must not be described as the +physical MITC4 strain field. + +### 19.2 FEM wiki cross-checks + +- `[[MITC4 Shell Element]]`: four-node degenerated-continuum kinematics, two + director-tangent rotations, and edge-midpoint assumed shear. +- `[[MITC Shell Kinematics]]`: midsurface/director/thickness interpolation. +- `[[Assumed Transverse Shear Strain Interpolation]]`: locking mechanism and tying + rationale. +- `[[Continuum Mechanics Based Four-Node Shell Element]]`: three-dimensional virtual + work lineage and benchmark thread. +- `[[Shell Locking Phenomenon]]`: separate transverse-shear and curved-shell membrane + locking risks. +- `[[Total Lagrangian Shell Formulation]]` and + `[[Green-Lagrange Strain Linearization]]`: future residual and material/geometric + tangent separation. + +Wiki pages are navigation and cross-check evidence; the approved requirements and +research brief remain the project source of truth. + +## 20. Requirement traceability + +| 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 | +| `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 | +| `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 | +| `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 | + +## 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. +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. + +Until these items are approved, the formulation is ready for Numerical Review but +not for Implementation Planning. + +### 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. +- 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. +- Preserve local frame and natural-coordinate identity without location averaging. + +### 21.3 Reference Model handoff + +- Use at least one S4 and one S4R source artifact, but compare formulation-independent + global displacement evidence rather than claiming element equivalence. +- Prioritize pinched cylinder and NAFEMS LE3 models compatible with nodal loads and + approved BC semantics. +- Propose the mixed U tolerance and nonblocking UR warning threshold. +- Do not create, repair, rename, or run reference artifacts during this formulation + gate. + +### 21.4 Implementation Planning handoff + +- Do not begin until Numerical Review closes Section 21.1. +- Translate the deterministic preprocessing, 24-to-20 transform, tying projection, + physical/drill split, quadrature, recovery, and invariant portfolio into + `RED -> GREEN -> VERIFY` tests before production changes. +- Keep future nonlinear state and tangent out of the current linear-static plan. + +### 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.