FESADev/docs/nonlinear-static-newton-raphson/formulation.md

Nonlinear Static Newton–Raphson 정식화

1. Metadata

• feature_id: nonlinear-static-newton-raphson
• source_requirement: docs/nonlinear-static-newton-raphson/requirements.md
• source_requirement_status: approved
• source_requirement_approval_date: 2026-08-23
• source_requirement_sha256: 910E351F3F90FFA42F07D15206B55E62C05EAE821E250424430148337506EBF2
• source_research: docs/nonlinear-static-newton-raphson/research.md
• source_research_status: passed
• source_research_sha256: C920A6B02381230D1DB419D73FE767D611B36720D27D6410507BA88C3CA28EEE
• status: ready-for-numerical-review
• owner_agent: formulation-agent
• date: 2026-08-23
• formulation_scope: fixed-initial-reference Total Lagrangian, materially linear and geometrically nonlinear original-style FESA-MITC4, Full Newton
• independence: FESA는 독립 solver이며 Abaqus, Midas 또는 다른 solver의 요소, 회전, stabilization, 적분 또는 상태 갱신 동등성을 주장하지 않는다.

이 문서는 Stage 3의 수학·알고리즘 계약이다. C++ API, parser, HDF5 path/schema, exact reference artifact, benchmark blocking 분류 및 수치 검토 tolerance는 정의하지 않는다. 기존 docs/linear-static-mitc4-shell/formulation.md의 §1–14에서 승인된 node order, initial director, local sign, MITC tying, constitutive boundary, 2 x 2 x 2 integration과 physical recovery를 보존한다. 그 문서 §15는 gap inventory와 conditional identity로만 사용하며 이 문서의 nonlinear 식을 대신하지 않는다.

2. Scope and Assumptions

2.1 Analysis boundary

• 해석 종류는 정확히 하나의 nonlinear static step이다.
• 하중 제어는 dimensionless scalar factor (\lambda\in[0,1])의 fixed increment path다.
• 각 증분은 Full Newton으로 평형을 찾는다. Modified Newton, tangent reuse, line search, damping, cutback, retry 및 arc-length는 없다.
• 모든 residual, tangent, load-work derivative와 convergence quantity는 같은 current trial identity에서 평가한다.
• formulation은 한 고정 initial reference configuration (\Omega_0)만 사용한다. accepted increment 뒤에도 reference geometry를 재설정하지 않는다.
• inertia, mass, damping 및 time integration은 N/A다.

2.2 Element and material boundary

• Source S4와 S4R은 같은 four-node original-style FESA-MITC4를 선택한다.
• Element는 Reissner–Mindlin-type degenerated continuum shell이며, thickness stretch와 (E_{33}), (S_{33})은 제외한다.
• Material은 centered homogeneous isotropic linear elasticity다.
• Geometry는 large displacement와 finite rotation을 허용하지만 material applicability는 small elastic strain이다. Material history/state variable은 없다.
• MITC4+ membrane field, reduced integration, hourglass control, EAS, physical drilling constitutive law, follower load 및 distributed load는 포함하지 않는다.

2.3 Coordinates, signs and units

• Natural coordinates are ((\xi,\eta,\zeta)\in[-1,1]^3).
• Node layout is (1=(-1,-1)), (2=(+1,-1)), (3=(+1,+1)), (4=(-1,+1)).
• (+\zeta) and the initial director define the positive shell face.
• Global nodal DOF order is ([UX,UY,UZ,URX,URY,URZ]).
• Translation and nodal force are global and have dimensions length and force.
• Rotation coordinates/corrections are dimensionless. Nodal moments have dimension force-length.
• Constitutive strain, second Piola–Kirchhoff stress and physical shell recovery use the fixed initial local frame at each reference midsurface location. No rotating output frame is inserted into the material tangent.
• User input uses one consistent unit system. FESA does not infer or convert units.

3. Primary Variables, DOFs and State

3.1 Fixed reference data

At node (I), let

[ \mathbf X_I\in\mathbb R^3,\qquad \mathbf Q_I^0= \begin{bmatrix}\mathbf a_I^0&\mathbf b_I^0&\mathbf d_I^0\end{bmatrix} \in SO(3) ]

be the initial midsurface position and the approved initial right-handed triad. Initial directors and tangent frames are constructed exactly by the approved linear MITC4 deterministic incident-element rules. Thickness (t>0), (E>0), and (-1<\nu<0.5) are fixed.

3.2 Global nonlinear coordinates

The nodal nonlinear coordinate is

[ \mathbf q_I= \begin{bmatrix}\mathbf u_I\\boldsymbol\phi_I\end{bmatrix} \in\mathbb R^6. ]

(\mathbf u_I) is global midsurface displacement. (\boldsymbol\phi_I) is a fixed-initial, fixed-global-basis rotation vector, and

[ \mathbf R_I(\boldsymbol\phi_I) =\exp([\boldsymbol\phi_I]_\times),\qquad \mathbf Q_I=\mathbf R_I\mathbf Q_I^0,\qquad \mathbf d_I=\mathbf Q_I\mathbf e_3 . ]

The mathematical branch is the principal exponential chart (|\boldsymbol\phi_I|<\pi). A target or trial exactly on the non-unique (\pi)-rotation branch, an ambiguous logarithm, or nonfinite rotation data fails closed. Numerical branch-distance tolerance is owned by Numerical Review; this formulation introduces no hidden clamp.

The 24-vector for one element is

[ \mathbf q_g= \begin{bmatrix} \mathbf u_1&\boldsymbol\phi_1& \cdots& \mathbf u_4&\boldsymbol\phi_4 \end{bmatrix}^{T}\in\mathbb R^{24}. ]

Constrained and free identities are the fixed global components of this vector.

3.3 Multiplicative update and retraction

For an unchanged Newton correction (\Delta\boldsymbol\phi_I), define

[ \mathbf Q_I^+ =\exp([\boldsymbol\phi_I+\Delta\boldsymbol\phi_I]_\times)\mathbf Q_I^0 =\Delta\mathbf R_I,\mathbf Q_I, ]

[ \boxed{ \Delta\mathbf R_I= \exp([\boldsymbol\phi_I+\Delta\boldsymbol\phi_I]\times) \exp(-[\boldsymbol\phi_I]\times)\in SO(3).} ]

Thus the update is left multiplication in the spatial/global frame. The correction is dimensionless, expressed in the fixed global basis, and is applied without a scalar modification. The authoritative orientation is (\mathbf Q_I); after every update (\boldsymbol\phi_I), (\mathbf Q_I), and (\mathbf d_I) must satisfy the boxed identity. Corrections are not multiplied successively as independent incremental rotation tensors; the fixed-reference exponential coordinate is updated first and (\Delta\mathbf R_I) is then derived. Therefore a final (\boldsymbol\phi_I) defines one path-independent triad inside the principal chart.

Rodrigues' map is

[ \exp([\boldsymbol\phi]\times) =\mathbf I+\frac{\sin\theta}{\theta}[\boldsymbol\phi]\times +\frac{1-\cos\theta}{\theta^2}[\boldsymbol\phi]_\times^2, \qquad \theta=|\boldsymbol\phi|, ]

with its analytic series at (\theta=0). No normalization of (\mathbf Q) is part of the mathematical update. A finite implementation may apply a checked orthogonal projection only as a representation safeguard if Numerical Review approves its invariant tolerance; it must not change residual/tangent state identity.

3.4 Exact first and second rotation variations

For a direction (\mathbf p), define the left Jacobian

[ \mathbf J_l(\boldsymbol\phi)= \mathbf I+ \frac{1-\cos\theta}{\theta^2}[\boldsymbol\phi]\times+ \frac{\theta-\sin\theta}{\theta^3}[\boldsymbol\phi]\times^2 . ]

Then the spatial infinitesimal rotation induced by the additive global-coordinate variation (\mathbf p) is

[ \boldsymbol\omega_{\mathbf p} =\mathbf J_l(\boldsymbol\phi)\mathbf p, \qquad D\mathbf d[\mathbf p] =\boldsymbol\omega_{\mathbf p}\times\mathbf d . ]

For directions (\mathbf p,\mathbf s),

[ \boxed{ D^2\mathbf d[\mathbf p,\mathbf s] =D\mathbf J_l(\boldsymbol\phi)[\mathbf s]\mathbf p\times\mathbf d +\boldsymbol\omega_{\mathbf p}\times (\boldsymbol\omega_{\mathbf s}\times\mathbf d).} ]

The symmetric value is obtained by the analytically identical exchange (\mathbf p\leftrightarrow\mathbf s); disagreement is an implementation error. (D\mathbf J_l) is the direct derivative of the displayed analytic expression, including scalar coefficient derivatives and the series branch. It is not omitted.

Equivalently, all required derivatives may be evaluated through the exact Fréchet primitives

[ D\exp_{\mathbf A}[\mathbf H] =\int_0^1 \exp((1-s)\mathbf A)\mathbf H\exp(s\mathbf A),ds, ]

[ D^2\exp_{\mathbf A}[\mathbf H,\mathbf K] =D!\left(D\exp_{\mathbf A}[\mathbf H]\right)[\mathbf K], ]

where the second expression is evaluated symmetrically. These definitions, the Rodrigues form, and the directional-derivative verification must give the same first/second variation.

3.5 Trial and committed state

A complete state identity contains at least

[ \mathcal S= {\lambda,\mathbf u,\boldsymbol\phi,\mathbf Q,\mathbf d, \mathbf f_{int},\mathbf f_{ext},\mathbf r, \text{reaction},\text{physical recovery rows}}. ]

Committed and trial instances never alias. Element drilling quantities are derived from the same trial and introduce no material history. A successful increment atomically transfers the entire validated trial identity. Nonconvergence, nonfinite data, invalid geometry/rotation, solve failure, recovery failure, or output failure discards the candidate and preserves the preceding committed state.

4. Strong Form and Boundary Conditions

4.1 Total Lagrangian parent continuum

On the fixed reference volume (\Omega_0),

[ \operatorname{Div}_0\mathbf P+\mathbf B_0=\mathbf0,\qquad \mathbf P=\mathbf F\mathbf S , ]

[ \mathbf F=\frac{\partial\mathbf x}{\partial\mathbf X},\qquad \mathbf E=\frac12(\mathbf F^T\mathbf F-\mathbf I),\qquad \mathbf S=\mathbb C_{ps}:\overline{\mathbf E}. ]

The parent boundary is split into (\Gamma_u\cup\Gamma_t):

[ \mathbf x=\overline{\mathbf x}\quad\text{on }\Gamma_u,\qquad \mathbf P\mathbf N_0=\overline{\mathbf T}_0\quad\text{on }\Gamma_t . ]

For this executable subset (\mathbf B_0=\mathbf0) and distributed (\overline{\mathbf T}_0=\mathbf0). External work is supplied only by approved nodal CLOAD and nodal prescribed global coordinates. The continuum natural term documents the lineage but does not authorize pressure, edge traction, body force, or follower loading.

4.2 Translational prescribed path

For a constrained translation component (c), step-start value (u_c^s), factor-one absolute target (u_c^\star), and current factor (\lambda),

[ \boxed{u_c(\lambda)=u_c^s+\lambda(u_c^\star-u_c^s).} ]

The step has no prior history in this release, but the formula preserves the step-start identity explicitly.

4.3 Finite rotational prescribed path

For a constrained global rotation-vector component (c),

[ \boxed{ \phi_c(\lambda)=\phi_c^s+ \lambda(\phi_c^\star-\phi_c^s).} ]

(\boldsymbol\phi^\star) is the factor-one absolute global rotation-vector target, not a Newton correction. At an increment target, constrained coordinates are set to the boxed absolute values before the first residual evaluation. Every Newton solve then uses

[ \boxed{\Delta\phi_c=0.} ]

Thus target advance and Newton correction are distinct. The triad path is (\mathbf Q(\lambda)=\exp([\boldsymbol\phi(\lambda)]_\times)\mathbf Q^0) only for components fixed by boundary data; free components remain solution variables. This componentwise path is deterministic, reaches the exact source target at (\lambda=1), and introduces no follower or controlled-displacement semantics.

4.4 Nodal load path and drilling-moment rejection

Approved nodal force and immutable source-moment coefficients (\mathbf F_I^\star,\mathbf M_{I,src}^\star) remain fixed in global axes:

[ \mathbf F_I(\lambda)=\lambda\mathbf F_I^\star,\qquad \mathbf M_{I,src}(\lambda)=\lambda\mathbf M_{I,src}^\star . ]

The factor changes magnitude only. A zero aggregated nodal moment is accepted as a separate exact-zero case. Every nonzero moment is checked against the same current trial director used by the physical residual:

[ \boxed{ \rho_{M,I}= \frac{|\mathbf d_I\cdot\mathbf M_{I,src}^\star|} {|\mathbf M_{I,src}^\star|} \le 10^{-12}.} ]

Failure at preprocessing or at any entered trial is unsupported-drilling-load and occurs before that trial residual/tangent can be accepted. The denominator is never clamped. This admissibility check does not project, rotate, shorten, overwrite, or otherwise mutate the source vector; source identity and the proportional fixed-global path remain exact.

The shell has no physical director-parallel work channel. After admissibility, the exact covector representing the immutable source moment on the current physical two-DOF virtual-rotation space is

[ \boxed{ \mathbf M_{I,phys}(\mathbf q;\lambda) =\lambda\mathbf P_I(\mathbf q)\mathbf M_{I,src}^\star,\qquad \mathbf P_I=\mathbf I-\mathbf d_I\mathbf d_I^T.} ]

Thus an admitted source component within the approved (10^{-12}) direction band is removed from physical work exactly and is not sent to the numerical drill channel. (\mathbf M_{I,phys}) is an effective work covector, not a mutated source load. Its state dependence is solely the restriction to the changing admissible director-tangent space; it does not rotate the fixed-global source direction and does not introduce follower-load semantics. Residual and tangent must include (D\mathbf P_I) exactly as specified in §11.2.

5. Weak or Variational Form

Let (\mathcal V_{24}(\lambda)) be the admissible trial manifold satisfying the absolute constraints and let (\mathcal V_{24}^0) contain variations with zero constrained global components. The equilibrium virtual-work statement is

[ \boxed{ \delta\mathcal G(\mathbf q_g;\lambda) =\delta\Pi_{phys} +\delta\Pi_{drill} -\delta W_F -\delta W_M=0 \quad\forall\delta\mathbf q_g\in\mathcal V_{24}^0 .} ]

The physical internal virtual work is

[ \delta\Pi_{phys} =\int_{\Omega_0} \delta\overline{\mathbf e}^{T}\mathbf s,dV_0, \qquad \mathbf s=\mathbf C_5\overline{\mathbf e}. ]

The translational dead-load work is

[ W_F=\lambda\sum_I\mathbf F_I^\star\cdot\mathbf u_I . ]

For a virtual rotation-coordinate variation (\delta\boldsymbol\phi_I), define the induced spatial rotation and its physical director-tangent part by

[ \delta\boldsymbol\vartheta_I =\mathbf J_l(\boldsymbol\phi_I)\delta\boldsymbol\phi_I,\qquad \delta\boldsymbol\vartheta_I^{phys} =(\mathbf I-\mathbf d_I\mathbf d_I^T) \delta\boldsymbol\vartheta_I =\mathbf d_I\times\delta\mathbf d_I . ]

[ \boxed{ \delta W_M =\sum_I\mathbf M_{I,phys}(\mathbf q;\lambda)\cdot \delta\boldsymbol\vartheta_I =\lambda\sum_I\mathbf M_{I,src}^\star\cdot \delta\boldsymbol\vartheta_I^{phys}.} ]

Here (\mathbf M_{I,src}^\star) is the immutable global source vector and (\mathbf M_{I,phys}=\lambda\mathbf P_I\mathbf M_{I,src}^\star) is the exact effective physical-work covector from §4.4. The two expressions are identical because (\delta\boldsymbol\vartheta_I^{phys}=\mathbf P_I \delta\boldsymbol\vartheta_I). Admissibility permits a source component within the approved ratio band; it does not assert (\mathbf M_{I,src}^\star\cdot\mathbf d_I=0) exactly. Projection removes that component from physical work, and the numerical director-parallel rotation does no external work. Configuration dependence of the generalized load arises from the exact rotation-coordinate and projector derivatives; the immutable source vector itself is never projected in storage, rotated, or shortened.

The physical shell and drilling terms are conservative inside their valid charts. An immutable constant-spatial source moment restricted to the physical tangent space as above defines a work one-form on finite rotations, not a global scalar rotation potential. Its exact load-work tangent is therefore generally nonsymmetric even though it contains no follower direction update. With no applied moment, the complete unreduced internal/drill tangent is a symmetric Hessian in exact arithmetic. Numerical symmetry expectations and norms belong to Numerical Review.

6. Discretization

6.1 Bilinear interpolation

[ N_I(\xi,\eta)= \frac14(1+\xi_I\xi)(1+\eta_I\eta), ]

[ \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} ]

Mandatory algebraic checks are

[ \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. ]

These checks apply independently of deformation.

6.2 Reference and current shell geometry

With (z=t\zeta/2),

[ \boxed{ \mathbf X(\xi,\eta,\zeta) =\sum_I N_I\mathbf X_I+ z\sum_I N_I\mathbf d_I^0,} ]

[ \boxed{ \mathbf x(\xi,\eta,\zeta;\mathbf q_g) =\sum_I N_I(\mathbf X_I+\mathbf u_I)+ z\sum_I N_I\mathbf d_I.} ]

Straight material lines through thickness remain straight. The current director is unit length by the rotation map. Thickness is constant and is not a state variable.

The continuum director field used by geometry is

[ \mathbf D(\xi,\eta)=\sum_I N_I\mathbf d_I . ]

It is bilinearly interpolated and is not normalized before entering (\mathbf x), (\mathbf g_\alpha), strain, tying, residual, or tangent. This is the selected original degenerated-continuum convention; silently replacing it by (\mathbf D/|\mathbf D|) changes the kinematics. Its exact variations are simply

[ D\mathbf D=\sum_I N_I D\mathbf d_I,\qquad D^2\mathbf D=\sum_I N_I D^2\mathbf d_I. ]

When a unit interpolated director is needed only for a checked diagnostic/current frame, use (\overline{\mathbf d}=\mathbf D/|\mathbf D|) after requiring (|\mathbf D|>0). With (\mathbf n=\mathbf v/|\mathbf v|) and (\mathbf P=\mathbf I-\mathbf n\mathbf n^T),

[ D\mathbf n[\mathbf p]=\frac{\mathbf P\mathbf p}{|\mathbf v|}, ]

[ \boxed{ D^2\mathbf n[\mathbf p,\mathbf s] =-\frac{ (\mathbf n\cdot\mathbf s)\mathbf P\mathbf p+ (\mathbf n\cdot\mathbf p)\mathbf P\mathbf s+ \mathbf n(\mathbf p\cdot\mathbf P\mathbf s)} {|\mathbf v|^2}.} ]

No normalized diagnostic frame feeds back into the physical constitutive law. At every trial require

[ \mathbf d_I\cdot\mathbf d_J>-1\quad(I<J), \qquad |\mathbf D(\xi,\eta)|>0 ]

at every required tying, Gauss, center and recovery location. Thus an exactly opposed nodal pair or a degenerate interpolated director fails closed rather than choosing a replacement axis. Scale-aware proximity tolerances, if required, belong to Numerical Review.

Define

[ \mathbf G_\alpha=\mathbf X_{,\alpha},\qquad \mathbf g_\alpha=\mathbf x_{,\alpha},\qquad \alpha\in{\xi,\eta,\zeta}. ]

Constant translation is reproduced exactly. For any proper rigid (\mathbf R) and translation (\mathbf c),

[ \mathbf x_I=\mathbf R\mathbf X_I+\mathbf c,\qquad \mathbf Q_I=\mathbf R\mathbf Q_I^0 ]

gives (\mathbf g_\alpha=\mathbf R\mathbf G_\alpha) and zero physical strain.

7. Kinematics and the Nonlinear 24-to-20 Map

7.1 Fixed-reference strain

[ \mathbf F=\mathbf g_\alpha\otimes\mathbf G^\alpha,\qquad \mathbf G^\alpha\cdot\mathbf G_\beta=\delta^\alpha_\beta, ]

[ \boxed{ E_{\alpha\beta} =\frac12(\mathbf g_\alpha\cdot\mathbf g_\beta -\mathbf G_\alpha\cdot\mathbf G_\beta).} ]

This is Green–Lagrange strain. Second Piola–Kirchhoff stress is its work-conjugate stress.

7.2 Exact physical director chart

Let

[ \mathbf c_I=\mathbf d_I^0\times\mathbf d_I,\qquad s_I=|\mathbf c_I|,\qquad c_I=\mathbf d_I^0\cdot\mathbf d_I. ]

For (c_I>-1), define the minimal swing vector

[ \boxed{ \boldsymbol\chi_I= \begin{cases} \mathbf0,&s_I=0,\ c_I=1,\[1mm] \dfrac{\operatorname{atan2}(s_I,c_I)}{s_I}\mathbf c_I, &s_I>0 . \end{cases}} ]

(\boldsymbol\chi_I\cdot\mathbf d_I^0=0) and (\exp([\boldsymbol\chi_I]_\times)\mathbf d_I^0=\mathbf d_I). The antipodal state (c_I=-1) is not a unique two-DOF director chart and fails closed.

Define the two physical rotation coordinates

[ \boldsymbol\chi_I^{2}= \begin{bmatrix} \mathbf a_I^0\cdot\boldsymbol\chi_I\ \mathbf b_I^0\cdot\boldsymbol\chi_I \end{bmatrix}\in\mathbb R^2 ]

and

[ \boxed{ \mathbf q_{20}=\boldsymbol\Phi(\mathbf q_g)

\begin{bmatrix} \mathbf u_1&\boldsymbol\chi_1^2& \cdots& \mathbf u_4&\boldsymbol\chi_4^2 \end{bmatrix}^{T}\in\mathbb R^{20}.} ]

This map discards only the rotation about the director gauge. It is an exact finite map, not the unresolved linear projection from the earlier gap inventory.

7.3 Map Jacobian and Hessian

Define

[ \mathbf A(\mathbf q_g) =D\boldsymbol\Phi(\mathbf q_g) \in\mathbb R^{20\times24}, \qquad \mathbf H_a^\Phi(\mathbf q_g) =D^2\Phi_a(\mathbf q_g) \in\mathbb R^{24\times24}. ]

The translation blocks are identity and have zero Hessian. Rotation blocks are obtained by the exact chain

[ \boldsymbol\phi \xrightarrow{\exp} \mathbf d \xrightarrow{\text{minimal swing}} \boldsymbol\chi \xrightarrow{(\mathbf a^0,\mathbf b^0)} \boldsymbol\chi^2 . ]

For implementation, first use §3.4 for (D\mathbf d,D^2\mathbf d), then differentiate (\mathbf c=\mathbf d^0\times\mathbf d), (s=|\mathbf c|), (c=\mathbf d^0\cdot\mathbf d), and (\alpha=\operatorname{atan2}(s,c)). The exact primitives are

[ D|\mathbf v|[\mathbf p] =\frac{\mathbf v\cdot\mathbf p}{|\mathbf v|}, ]

[ D^2|\mathbf v|[\mathbf p,\mathbf s] =\frac{\mathbf p\cdot\mathbf s}{|\mathbf v|} -\frac{(\mathbf v\cdot\mathbf p)(\mathbf v\cdot\mathbf s)} {|\mathbf v|^3}, ]

[ D\operatorname{atan2}(y,x)[p] =\frac{x,Dy[p]-y,Dx[p]}{x^2+y^2}, ]

with the direct derivative of this last expression for the mixed second variation. The analytic (s\rightarrow0,c\rightarrow1) series must be used; neither first nor second derivative may be obtained by normalizing a zero cross product.

At the undeformed state,

[ \boxed{ \mathbf A_0=\mathbf T_p,} ]

where (\mathbf T_p) is the approved linear 20-by-24 physical transform, and

[ D\boldsymbol\chi_I^2[\delta\boldsymbol\phi_I]

\begin{bmatrix} (\mathbf a_I^0)^T\ (\mathbf b_I^0)^T \end{bmatrix} \delta\boldsymbol\phi_I . ]

The matrices (\mathbf H_a^\Phi) are generally nonzero away from the linear state. They are required in the global tangent.

7.4 First and second strain variations

For a physical coordinate (q_{20,A}), define

[ \mathbf h_A=\frac{\partial\mathbf x}{\partial q_{20,A}}, \qquad \mathbf h_{AB}= \frac{\partial^2\mathbf x} {\partial q_{20,A}\partial q_{20,B}}. ]

Then

[ \boxed{ E_{\alpha\beta,A} =\frac12( \mathbf h_{A,\alpha}\cdot\mathbf g_\beta+ \mathbf g_\alpha\cdot\mathbf h_{A,\beta}),} ]

[ \boxed{ \begin{aligned} E_{\alpha\beta,AB} =\frac12(& \mathbf h_{AB,\alpha}\cdot\mathbf g_\beta+ \mathbf g_\alpha\cdot\mathbf h_{AB,\beta}\ &+\mathbf h_{A,\alpha}\cdot\mathbf h_{B,\beta} +\mathbf h_{B,\alpha}\cdot\mathbf h_{A,\beta}). \end{aligned}} ]

Director-coordinate blocks retain the exact first and second derivative of (\exp([\boldsymbol\chi]\times)\mathbf d^0). Translation-translation (\mathbf h{AB}) is zero, but the quadratic gradient products remain and generate stress-dependent geometric stiffness.

8. MITC4 Tying and Local Engineering Measures

8.1 Tying locations

[ T_{\xi-}=(0,-1,0),\quad T_{\xi+}=(0,+1,0),\quad T_{\eta-}=(-1,0,0),\quad T_{\eta+}=(+1,0,0). ]

The first pair supplies (E_{\xi\zeta}) and the second pair supplies (E_{\eta\zeta}).

8.2 Nonlinear assumed shear

[ \boxed{ \overline E_{\xi\zeta}(\xi,\eta) =\frac{1-\eta}{2}E_{\xi\zeta}^{DI}(T_{\xi-}) +\frac{1+\eta}{2}E_{\xi\zeta}^{DI}(T_{\xi+}),} ]

[ \boxed{ \overline E_{\eta\zeta}(\xi,\eta) =\frac{1-\xi}{2}E_{\eta\zeta}^{DI}(T_{\eta-}) +\frac{1+\xi}{2}E_{\eta\zeta}^{DI}(T_{\eta+}).} ]

All in-plane components remain direct. The same linear projection operator is applied to the strain and both variations:

[ \overline{\mathbf E}=\mathcal P_{MITC}\mathbf E,\qquad D\overline{\mathbf E}=\mathcal P_{MITC}D\mathbf E,\qquad D^2\overline{\mathbf E}=\mathcal P_{MITC}D^2\mathbf E . ]

Tying geometry/directors and the trial identity must be identical for residual and tangent. Direct shear in any one of the three is prohibited.

8.3 Fixed reference local frame

At each reference midsurface location, construct the approved initial frame ((\mathbf e_1^0,\mathbf e_2^0,\mathbf e_3^0)). Transform the assumed tensor into that frame and define

[ \boxed{ \overline{\mathbf e}= \begin{bmatrix} E_{11}&E_{22}&2E_{12}&2E_{13}&2E_{23} \end{bmatrix}^{T}.} ]

Let

[ \mathbf B_A=\frac{\partial\overline{\mathbf e}} {\partial q_{20,A}}\in\mathbb R^5,\qquad \mathbf G_{AB}= \frac{\partial^2\overline{\mathbf e}} {\partial q_{20,A}\partial q_{20,B}}\in\mathbb R^5 . ]

The fixed frame has no trial derivative. A rotating recovery frame, if ever defined downstream, cannot be substituted into this constitutive calculation.

9. Constitutive Contract

Let

[ G=\frac{E}{2(1+\nu)},\qquad \kappa_s=\frac56, ]

[ \mathbf C_{ps}= \frac{E}{1-\nu^2} \begin{bmatrix} 1&\nu&0\ \nu&1&0\ 0&0&(1-\nu)/2 \end{bmatrix}, ]

[ \boxed{ \mathbf C_5= \begin{bmatrix} \mathbf C_{ps}&\mathbf0\ \mathbf0&\kappa_sG\mathbf I_2 \end{bmatrix}.} ]

The second Piola–Kirchhoff engineering stress vector is

[ \boxed{ \mathbf s= \begin{bmatrix} S_{11}&S_{22}&S_{12}&S_{13}&S_{23} \end{bmatrix}^{T} =\mathbf C_5\overline{\mathbf e}.} ]

(\mathbf C_T=\mathbf C_5) is constant, symmetric and positive definite for the approved material range. This is a St. Venant–Kirchhoff-type shell restriction: material response is linear in Green–Lagrange strain and is not evidence for large material strain validity. There are no constitutive updates or state variables.

10. Objective Numerical Drilling Contract

10.1 Element corotation

At the element center, form

[ \mathbf J_0^c= \begin{bmatrix}\mathbf G_\xi&\mathbf G_\eta&\mathbf G_\zeta\end{bmatrix}{(0,0,0)}, \qquad \mathbf j^c= \begin{bmatrix}\mathbf g\xi&\mathbf g_\eta&\mathbf g_\zeta\end{bmatrix}_{(0,0,0)}, ]

[ \mathbf F_c=\mathbf j^c(\mathbf J_0^c)^{-1}. ]

Require finite (\det\mathbf F_c>0). Its right polar decomposition is

[ \mathbf F_c=\mathbf R_e\mathbf U_e,\qquad \mathbf R_e\in SO(3),\quad \mathbf U_e=(\mathbf F_c^T\mathbf F_c)^{1/2}\succ0 . ]

For a superposed rigid motion, (\mathbf F_c^+=\mathbf R\mathbf F_c) and (\mathbf R_e^+=\mathbf R\mathbf R_e).

10.2 Objective drill angle

The nodal total triad rotation is

[ \mathbf R_I=\mathbf Q_I(\mathbf Q_I^0)^T . ]

Remove the element corotation:

[ \mathbf C_{eI}=\mathbf R_e^T\mathbf R_I. ]

Let (\mathbf y_{eI}=\mathbf C_{eI}\mathbf d_I^0). Construct the minimal swing (\mathbf S_{eI}) that maps (\mathbf d_I^0) to (\mathbf y_{eI}) using the same swing formula as §7.2, and define

[ \mathbf T_{eI}=\mathbf S_{eI}^T\mathbf C_{eI}. ]

(\mathbf T_{eI}\mathbf d_I^0=\mathbf d_I^0), so it is a twist about the initial director. With (\operatorname{vex}([\mathbf v]_\times)=\mathbf v),

[ s_{\gamma,I}= \frac12\mathbf d_I^0\cdot \operatorname{vex}(\mathbf T_{eI}-\mathbf T_{eI}^T), \qquad c_{\gamma,I}=\frac{\operatorname{tr}\mathbf T_{eI}-1}{2}, ]

[ \boxed{ \gamma_{eI}=\operatorname{atan2}(s_{\gamma,I},c_{\gamma,I}).} ]

The vector (\boldsymbol\gamma_e=[\gamma_{e1},\ldots,\gamma_{e4}]^T\in\mathbb R^4) contains numerical coordinates only. Antipodal swing or the non-unique twist branch at (\pm\pi) fails closed.

For (\mathbf x=\mathbf R\mathbf X+\mathbf c) and (\mathbf Q_I=\mathbf R\mathbf Q_I^0), (\mathbf R_e=\mathbf R), (\mathbf C_{eI}=\mathbf I), and (\boldsymbol\gamma_e=\mathbf0). Therefore finite rigid motion creates no drilling energy.

10.3 Fixed scale and potential

Compute the approved linear physical (\mathbf K_{20}^0) once from the initial reference element. Let (\mathcal R_+) contain only its eight finite positive director-tangent rotational diagonals. Require (\mathcal R_+\ne\varnothing), then fix

[ \boxed{ k_{ref,e}=\min\mathcal R_+,\qquad k_{d,e}=10^{-3}k_{ref,e}.} ]

(k_{d,e}) has force-length units and remains immutable for all increments and iterations. It is not recomputed from a trial or committed tangent, so its residual and Hessian have no hidden scale derivative.

[ \boxed{ \Pi_{drill}^e=\frac12 k_{d,e} \boldsymbol\gamma_e^T\boldsymbol\gamma_e.} ]

It is not physical shell strain energy and creates no physical strain, stress, resultant, or drilling output.

10.4 Exact drill residual and Hessian

Define

[ \mathbf D_\gamma= D\boldsymbol\gamma_e\in\mathbb R^{4\times24}, \qquad \mathbf H_I^\gamma=D^2\gamma_{eI}\in\mathbb R^{24\times24}. ]

Then

[ \boxed{ \mathbf r_{drill}^e =k_{d,e}\mathbf D_\gamma^T\boldsymbol\gamma_e,} ]

[ \boxed{ \mathbf K_{drill}^e =k_{d,e}\left( \mathbf D_\gamma^T\mathbf D_\gamma+ \sum_{I=1}^4\gamma_{eI}\mathbf H_I^\gamma\right).} ]

Both derivatives include (\mathbf R_e), nodal exponential maps, swing removal and atan2. A frozen-frame or secant drill tangent is prohibited.

For exact polar derivatives, write (\delta\mathbf F=\mathbf R_e\mathbf A), (\delta\mathbf R_e=\mathbf R_e\boldsymbol\Omega), where (\boldsymbol\Omega^T=-\boldsymbol\Omega). Solve

[ \boxed{ \boldsymbol\Omega\mathbf U_e+\mathbf U_e\boldsymbol\Omega =\mathbf A-\mathbf A^T,} ]

and set (\delta\mathbf U_e=\mathbf A-\boldsymbol\Omega\mathbf U_e). For mixed directions (p,s), define

[ \mathbf A_p=\mathbf R_e^T\mathbf F_{c,p},\qquad \mathbf A_{p,s}=-\boldsymbol\Omega_s\mathbf A_p+ \mathbf R_e^T\mathbf F_{c,ps}, ]

and differentiate this Sylvester equation:

[ \boldsymbol\Omega_{ps}\mathbf U_e+ \mathbf U_e\boldsymbol\Omega_{ps} =\mathbf A_{p,s}-\mathbf A_{p,s}^T -\boldsymbol\Omega_p\mathbf U_{e,s} -\mathbf U_{e,s}\boldsymbol\Omega_p , ]

[ \mathbf R_{e,ps} =\mathbf R_e(\boldsymbol\Omega_s\boldsymbol\Omega_p+ \boldsymbol\Omega_{ps}). ]

The exchange (p\leftrightarrow s) must agree. Applying the norm, swing and atan2 derivative primitives from §7.3 completes exact (\mathbf D_\gamma,\mathbf H_I^\gamma).

At the initial state, a pure nodal drill direction (\delta\boldsymbol\phi_I=\gamma_I\mathbf d_I^0), with zero translations and zero physical swing, gives the approved (k_{d,e}\mathbf I_4) action. General deformation also contains the objective corotation coupling required to preserve the six physical rigid modes.

11. Element Residual and Same-Trial Tangent

11.1 Physical 20-coordinate equations

[ \Pi_{phys}^e =\frac12\int_{\Omega_0^e} \overline{\mathbf e}^T\mathbf C_5\overline{\mathbf e},dV_0. ]

The physical residual and tangent are

[ \boxed{ r_{20,A}^{phys} =\int_{\Omega_0^e}\mathbf B_A^T\mathbf s,dV_0,} ]

[ \boxed{ K_{20,AB}^{mat} =\int_{\Omega_0^e} \mathbf B_A^T\mathbf C_5\mathbf B_B,dV_0,} ]

[ \boxed{ K_{20,AB}^{geo} =\int_{\Omega_0^e} \mathbf s^T\mathbf G_{AB},dV_0,} ]

[ \mathbf r_{20}^{phys}\in\mathbb R^{20},\qquad \mathbf K_{20}^{phys} =\mathbf K_{20}^{mat}+\mathbf K_{20}^{geo} \in\mathbb R^{20\times20}. ]

The geometric term includes the exact director Hessian and the quadratic Green–Lagrange terms. It is not a separately guessed initial-stress matrix.

11.2 Fixed-spatial nodal load in physical coordinates

In the element equations below, (\mathbf F_I^\star,\mathbf M_{I,src}^\star) denote a deterministic source-ordered element load share whose incident-element sum is the single aggregated nodal CLOAD. Equivalently the identical nodal contribution may be assembled once outside element kernels. It must never be copied in full to every incident element. Define the per-unit-load-factor effective physical covector

[ \mathbf M_{I,phys}^\star(\mathbf q) :=\mathbf P_I(\mathbf q)\mathbf M_{I,src}^\star, \qquad \mathbf M_{I,phys}(\mathbf q;\lambda) =\lambda\mathbf M_{I,phys}^\star(\mathbf q). ]

For nodal physical swing coordinates (\chi_{I,a}^2,\ a\in{1,2}), define

[ \mathbf L_{I,a} =\mathbf d_I\times \frac{\partial\mathbf d_I}{\partial\chi_{I,a}^2} \in\mathbb R^3, \qquad m_{I,a}^{20} =\mathbf M_{I,phys}^\star\cdot\mathbf L_{I,a} =\mathbf M_{I,src}^\star\cdot\mathbf L_{I,a}. ]

Insert these two entries per node, and zeros in translation positions, into (\mathbf m_{20}^{M}(\mathbf q_{20})\in\mathbb R^{20}). Since (\mathbf d_I\cdot\mathbf L_{I,a}=0), the two displayed forms of (m_{I,a}^{20}) are exactly identical even when the admitted source has a small director component. Its exact physical-chart derivative is

[ \begin{aligned} K_{I,ab}^{M,20} &=\frac{\partial\mathbf M_{I,phys}^\star} {\partial\chi_{I,b}^2}\cdot\mathbf L_{I,a} +\mathbf M_{I,phys}^\star\cdot \frac{\partial\mathbf L_{I,a}}{\partial\chi_{I,b}^2}\ &=\mathbf M_{I,src}^\star\cdot\left( \frac{\partial\mathbf d_I}{\partial\chi_{I,b}^2} \times \frac{\partial\mathbf d_I}{\partial\chi_{I,a}^2} + \mathbf d_I\times \frac{\partial^2\mathbf d_I} {\partial\chi_{I,a}^2\partial\chi_{I,b}^2} \right). \end{aligned} ]

The first line explicitly contains (\partial\mathbf M_{I,phys}^\star/\partial\chi_{I,b}^2 =(\partial\mathbf P_I/\partial\chi_{I,b}^2) \mathbf M_{I,src}^\star). The second line is the algebraically identical immutable-source form obtained by differentiating (\mathbf M_{I,src}^\star\cdot\mathbf L_{I,a}); therefore the projector derivative must not be dropped when the first form is used.

All cross-node and translation derivatives are zero. In general (\mathbf K_{20}^{M}) need not be symmetric; this is the exact derivative of a constant-spatial source-moment work one-form restricted to physical swing space.

Let (\mathbf f_{24}^{F,\star}\in\mathbb R^{24}) be the translational dead-force selection vector. The element external generalized vector and its exact derivative are

[ \boxed{ \mathbf f_{ext}^e =\lambda\left( \mathbf f_{24}^{F,\star} +\mathbf A^T\mathbf m_{20}^{M}\right),} ]

[ \boxed{ \mathbf K_{ext}^{M,e} =\lambda\left[ \mathbf A^T\mathbf K_{20}^{M}\mathbf A +\sum_{a=1}^{20}m_{20,a}^{M}\mathbf H_a^\Phi \right].} ]

As an independent exact global-coordinate identity, the per-unit-load-factor nodal rotational block is

[ \mathbf f_{\phi,I}^{M,\star} =\mathbf J_l(\boldsymbol\phi_I)^T \mathbf P_I\mathbf M_{I,src}^\star =\mathbf J_l(\boldsymbol\phi_I)^T\mathbf M_{I,phys}^\star, \qquad \mathbf P_I=\mathbf I-\mathbf d_I\mathbf d_I^T, ]

[ D\mathbf f_{\phi,I}^{M,\star}[\mathbf p] =D\mathbf J_l(\boldsymbol\phi_I)[\mathbf p]^T \mathbf P_I\mathbf M_{I,src}^\star +\mathbf J_l(\boldsymbol\phi_I)^T D\mathbf P_I[\mathbf p]\mathbf M_{I,src}^\star, ]

[ D\mathbf P_I[\mathbf p] =-\left( D\mathbf d_I[\mathbf p]\mathbf d_I^T+ \mathbf d_I D\mathbf d_I[\mathbf p]^T\right). ]

Let (\mathbf f_{24}^{M,\star}) contain these four nodal rotational blocks and zero translation blocks. Virtual-work covariance of the exact map gives

[ \boxed{ \mathbf A^T\mathbf m_{20}^{M}=\mathbf f_{24}^{M,\star}, \qquad D(\mathbf A^T\mathbf m_{20}^{M})[\mathbf p] =\left(\mathbf A^T\mathbf K_{20}^{M}\mathbf A +\sum_{a=1}^{20}m_{20,a}^{M}\mathbf H_a^\Phi\right)\mathbf p =D\mathbf f_{24}^{M,\star}[\mathbf p].} ]

Thus the physical-chart and global-coordinate forms are the same residual covector and the same derivative of the §5 work one-form. Both include the projector derivative; neither changes the immutable source identity nor creates follower-load direction updates or a drilling load channel.

11.3 Complete 24-coordinate equations

The element residual is

[ \boxed{ \mathbf r_e= \mathbf A^T\mathbf r_{20}^{phys} +\mathbf r_{drill}^e -\lambda\left( \mathbf f_{24}^{F,\star} +\mathbf A^T\mathbf m_{20}^{M}\right) \in\mathbb R^{24}.} ]

The exact same-trial algorithmic tangent is

[ \boxed{ \begin{aligned} \mathbf K_{T,e} =&\ \mathbf A^T\mathbf K_{20}^{phys}\mathbf A\ &+\sum_{a=1}^{20}r_{20,a}^{phys}\mathbf H_a^\Phi +\mathbf K_{drill}^e\ &-\lambda\left[ \mathbf A^T\mathbf K_{20}^{M}\mathbf A +\sum_{a=1}^{20}m_{20,a}^{M}\mathbf H_a^\Phi \right] \in\mathbb R^{24\times24}. \end{aligned}} ]

The first line contains mapped material and Green–Lagrange geometric terms. The second line contains physical coordinate-map curvature and the objective drill gradient/Hessian. The last two lines are the complete fixed-spatial source-moment work-conjugacy derivative, including the effective-work projector derivative, rotation-coordinate derivative and map-Hessian term. Fixed-global nodal force has no configuration derivative.

No term may be evaluated from a committed or previous-iteration state while another term uses the current trial. In particular, the following same-trial identities are mandatory:

[ \mathbf K_{T,e}=D\mathbf r_e,\qquad \mathbf K_T=D\mathbf r,\qquad \mathbf r=\sum_e\mathcal A_e^T\mathbf r_e . ]

11.4 Global equations and reactions

Assemble element contributions and nodal load shares in stable source order:

[ \mathbf r(\mathbf q;\lambda) =\mathbf f_{int}(\mathbf q) -\mathbf f_{ext}(\mathbf q;\lambda), \qquad \mathbf K_T=D_{\mathbf q}\mathbf r. ]

At a trial that satisfies its prescribed coordinates, Full Newton on free coordinates is

[ \boxed{ \mathbf K_{ff}^{(k)}\Delta\mathbf q_f^{(k)} =-\mathbf r_f^{(k)},\qquad \Delta\mathbf q_c^{(k)}=\mathbf0.} ]

The tangent is freshly factorized and solved once for every entered iteration. A valid (0\times0) free system is a trivial solve, not singular.

The full residual is preserved. Its constrained entries are generalized reactions conjugate to global ((\mathbf u,\boldsymbol\phi)); free entries are equilibrium evidence. Translation entries are spatial forces. For balance or global physical moment reporting, convert a rotational generalized covector by

[ \boxed{ \mathbf m_{sp,I} =\mathbf J_l(\boldsymbol\phi_I)^{-T}\mathbf r_{\phi,I}.} ]

This conversion is applied consistently to internal, external and reaction covectors; it does not alter the Newton equations.

12. Mapping and Numerical Integration

12.1 Reference Jacobian and derivative transformation

[ \mathbf J_0= \begin{bmatrix}\mathbf G_\xi&\mathbf G_\eta&\mathbf G_\zeta\end{bmatrix}, \qquad J_0=\det\mathbf J_0, ]

[ dV_0=J_0,d\xi,d\eta,d\zeta,\qquad \nabla_0(\cdot)=(\cdot)_{,\alpha}\otimes\mathbf G^\alpha . ]

Every derivative transformation uses the inverse/contravariant basis of this fixed reference Jacobian. The current mapping

[ \mathbf j= \begin{bmatrix}\mathbf g_\xi&\mathbf g_\eta&\mathbf g_\zeta\end{bmatrix}, \qquad j=\det\mathbf j ]

is not the integration measure but must remain finite and positive to exclude local inversion.

12.2 Required validation locations

Validate finite data, nonzero surface area, finite (J_0>0), and finite (j>0) at:

• all eight volume integration points;
• all four midsurface MITC tying points;
• the element center used by drilling corotation;
• every bottom/middle/top recovery point before commit.

Also require finite reciprocal bases, nonzero interpolated director, valid fixed local frame, (\det\mathbf F_c>0), SPD polar stretch, valid swing/twist branches, and (Q^TQ=I,\det Q=+1). A failed point is not omitted, averaged, clamped, or replaced.

12.3 Quadrature

Use tensor-product Gauss points

[ \xi,\eta,\zeta\in \left{-\frac1{\sqrt3},+\frac1{\sqrt3}\right},\qquad w_\xi=w_\eta=w_\zeta=1. ]

Thus each element uses four midsurface by two thickness points, eight volume evaluations. The following are integrated at all eight points with the same reference (J_0w):

• physical internal residual;
• material tangent;
• geometric tangent;
• physical strain-energy and physical recovery thickness moments.

The MITC tying values are evaluated at the four (\zeta=0) tying locations and used for both thickness points. The objective drill potential is an element-level four-coordinate numerical term and is not multiplied by volume quadrature.

S4 and S4R use this identical rule. One-point integration, reduced/selective integration, source-selected hourglass behavior, and higher-order production quadrature are prohibited.

13. Output Recovery

Only accepted committed states are recoverable authoritative frames. No trial Newton state is an authoritative result.

13.1 Nodal kinematics and reactions

At each source node:

• displacement ([U1,U2,U3]=\mathbf u_I), global, length;
• rotation ([UR1,UR2,UR3]=\boldsymbol\phi_I), fixed-initial global rotation-vector components, dimensionless;
• committed (\mathbf Q_I\in SO(3)) and (\mathbf d_I=\mathbf Q_I\mathbf e_3) as the kinematic state identity;
• translation reaction from constrained full residual entries, global force;
• rotational generalized reaction and its spatial conversion from §11.4, force-length.

For unambiguous recovery, let (\mathbf C_{u,I},\mathbf C_{\phi,I}\in\mathbb R^{3\times3}) be diagonal selectors with one on a prescribed global component and zero otherwise. With the accepted full-residual blocks (\mathbf r_{u,I},\mathbf r_{\phi,I}\in\mathbb R^3), define

[ \boxed{ \mathbf F_I^{react}:=\mathbf C_{u,I}\mathbf r_{u,I} \in\mathbb R^3, \qquad \mathbf m_{\phi,I}^{react}:=\mathbf C_{\phi,I}\mathbf r_{\phi,I} \in\mathbb R^3,} ]

[ \boxed{ \mathbf M_I^{react} :=\mathbf J_l(\boldsymbol\phi_I)^{-T}\mathbf m_{\phi,I}^{react} \in\mathbb R^3.} ]

(\mathbf F_I^{react}) has force units and is work-conjugate to (\delta\mathbf u_I). The generalized covector (\mathbf m_{\phi,I}^{react}) and spatial reaction moment (\mathbf M_I^{react}) have force-length units and satisfy (\mathbf m_{\phi,I}^{react}\cdot\delta\boldsymbol\phi_I =\mathbf M_I^{react}\cdot\delta\boldsymbol\vartheta_I). Unselected free residual components remain equilibrium evidence and are not relabelled as reactions. The symbol (\mathbf R_I) remains reserved exclusively for the rotation matrix defined in §§3.2 and 10.2.

Exact HDF5 projection is I/O Definition ownership. This section fixes the quantity, coordinate and unit meaning only.

13.2 Integration-point physical strain and stress

At each four midsurface Gauss locations and both thickness Gauss points recover (\overline{\mathbf e}) and (\mathbf s) from the exact accepted state, same MITC tying, fixed reference frame and material used by residual/tangent.

Generalized strains are thickness moments:

[ \boldsymbol\epsilon_0 =\frac1t\int_{-t/2}^{t/2} \begin{bmatrix}E_{11}&E_{22}&2E_{12}\end{bmatrix}^T dz, ]

[ \boldsymbol\kappa =\frac{12}{t^3}\int_{-t/2}^{t/2} z\begin{bmatrix}E_{11}&E_{22}&2E_{12}\end{bmatrix}^T dz, ]

[ \boldsymbol\gamma_0 =\frac1t\int_{-t/2}^{t/2} \begin{bmatrix}2E_{13}&2E_{23}\end{bmatrix}^T dz . ]

Their order is

[ [E11,E22,G12,K11,K22,K12,G13,G23]. ]

Membrane/shear strain is dimensionless and curvature has (1/\text{length}).

13.3 Section resultants and stresses

[ \mathbf N=\int_{-t/2}^{t/2} \begin{bmatrix}S_{11}&S_{22}&S_{12}\end{bmatrix}^T dz, ]

[ \mathbf M=\int_{-t/2}^{t/2} z\begin{bmatrix}S_{11}&S_{22}&S_{12}\end{bmatrix}^T dz, \qquad \mathbf Q_s=\int_{-t/2}^{t/2} \begin{bmatrix}S_{13}&S_{23}\end{bmatrix}^T dz . ]

The order remains

[ [N11,N22,N12,M11,M22,M12,Q13,Q23]. ]

(\mathbf N,\mathbf Q_s) have force/length units and (\mathbf M) has force units. Positive (z) and positive bending sign are the approved linear convention. At (z=-t/2,0,+t/2), recover only ([S11,S22,S12]) as direct evaluations. (S_{33}=0) is documented but not emitted.

There is no nodal extrapolation/averaging and no numerical drill stress, resultant, energy ratio, or output. Exact element, natural-coordinate, thickness-position, frame and source identity remain attached to every physical row.

13.4 Equilibrium identities

At an accepted state,

[ \mathbf r_f=\mathbf0 ]

in the exact equation; numerical acceptance tolerance is downstream. Define the applied physical load resultants at the accepted state by

[ \mathbf F_I^{appl}(\lambda):=\lambda\mathbf F_I^\star \in\mathbb R^3, \qquad \mathbf M_I^{appl}(\mathbf q;\lambda) :=\mathbf M_{I,phys}(\mathbf q;\lambda) =\lambda\mathbf P_I\mathbf M_{I,src}^\star \in\mathbb R^3. ]

Here (\mathbf F_I^{appl}) has force units and (\mathbf M_I^{appl}) has force-length units. For a fixed global origin (\mathbf X_o), physical resultant evidence uses current nodal positions (\mathbf x_I=\mathbf X_I+\mathbf u_I):

[ \boxed{ \mathbf e_F =\sum_I(\mathbf F_I^{appl}+\mathbf F_I^{react}) \in\mathbb R^3\ \text{[force]},} ]

[ \boxed{ \mathbf e_M=\sum_I[ (\mathbf x_I-\mathbf X_o)\times (\mathbf F_I^{appl}+\mathbf F_I^{react}) +\mathbf M_I^{appl}+\mathbf M_I^{react}] \in\mathbb R^3\ \text{[force-length]}.} ]

The balance uses the exact effective physical-work moment, while the immutable source (\mathbf M_{I,src}^\star) remains unchanged for source identity and the fixed-global load path. Its removed director component is neither a physical balance contribution nor a numerical drill load. The objective drilling potential contributes self-equilibrated nodal forces/moments and is included in full residual reactions; it has no separate physical output. Balance acceptance tolerance is owned by Numerical Review.

14. Newton and State Algorithm Pseudocode

preprocess in stable source order:
    validate source S4/S4R mapping, material, section, node order and initial geometry
    construct approved initial nodal triads Q0 and directors d0
    aggregate immutable fixed-global nodal CLOAD sources and validate initial moment rho_M
    validate all reference Gauss/tying/center/recovery mappings
    compute initial linear physical K20_0 for each element
    select k_ref from finite positive physical tangent-rotation diagonals
    freeze k_d = 1e-3 * k_ref for the entire analysis

initialize:
    committed = exact step-start state at factor 0
    intended factors follow the approved fixed increment and final shortening

for each intended factor lambda in deterministic order:
    trial = deep copy(committed)
    set every constrained translation and rotation coordinate to its absolute path value
    keep the preceding committed free coordinates as the initial free guess

    for iteration k = 1..100:
        derive every Q = Exp(phi^) * Q0 and d = Q * e3 from this trial
        validate finite state, principal charts, current geometry and source-moment rho_M

        for each element in stable source order:
            evaluate one trial identity at all Gauss/tying/center locations
            form MITC Green-Lagrange strain and exact first/second variations
            form r20_phys, K20_mat and K20_geo
            form exact Phi Jacobian A and all map Hessians H_Phi
            form polar corotation, objective gamma, D_gamma and H_gamma
            form r_drill and exact K_drill
            form fixed-force and projected physical-moment work plus exact P/J_l/A derivatives
            form complete r_e and K_T_e from §11.3

        deterministically assemble full r and K_T
        restrict to the stable free coordinate order
        freshly factorize K_ff for this iteration
        solve K_ff * Delta_q_f = -r_f; set Delta_q_c = 0
        reject any nonfinite residual, tangent, factorization or correction

        evaluate all four inclusive convergence gates from:
            current same-trial free residual, internal force and external force
            fresh same-iteration free correction
            current trial incremental solution relative to step-start committed state
        emit exactly one deterministic diagnostic summary for entered iteration k

        if all four gates pass:
            recover all physical rows from this unchanged current trial
            validate complete finite candidate and state/output identities
            atomically commit displacement, phi, Q, d, factor, forces, reaction and rows
            create exactly one accepted-frame candidate identity
            exit the iteration loop

        if k == 100:
            discard trial and fail without iteration 101, cutback or retry

        apply the complete fresh correction unchanged:
            u <- u + Delta_u
            phi <- phi + Delta_phi
            Q_new <- Exp((phi_new)^) * Q0 = Delta_R * Q_old
        continue with a new same-trial evaluation

on any increment/recovery/output failure:
    discard the complete trial/candidate
    preserve prior committed state and prior authoritative final result
    perform no replacement increment or hidden solver fallback

14.1 Four convergence channels

Let free translation and rotation-coordinate subvectors carry superscripts (u) and (\phi). Define the matching external residual-channel scales

[ S_{ext}^u=|\mathbf f_{ext,f}^u|\infty, \qquad S{ext}^\phi=|\mathbf f_{ext,f}^\phi|_\infty. ]

When the matching external scale is nonzero, residual channels retain

[ \eta_R^u= \frac{|\mathbf r_f^u|\infty} {\max(|\mathbf f{int,f}^u|\infty, |\mathbf f{ext,f}^u|_\infty)}, ]

[ \eta_R^\phi= \frac{|\mathbf r_f^\phi|\infty} {\max(|\mathbf f{int,f}^\phi|\infty, |\mathbf f{ext,f}^\phi|_\infty)} ]

for (S_{ext}^u>0) and for (S_{ext}^\phi>0), respectively. Each nonzero-external-scale channel passes iff its ratio is (\le10^{-2}). An exact-zero matching external scale instead selects the existing typed absolute residual gate:

[ \boxed{ S_{ext}^u=0 \ \Longrightarrow |\mathbf r_f^u|_\infty\le10^{-3} \quad\text{model-consistent force},} ]

[ \boxed{ S_{ext}^\phi=0 \ \Longrightarrow |\mathbf r_f^\phi|_\infty\le10^{-3} \quad\text{model-consistent force-length}.} ]

The exact-zero branch selector is the matching external scale itself, not (\max(|\mathbf f_{int,f}|\infty, |\mathbf f{ext,f}|_\infty)). Therefore a nonzero internal residual does not select the relative branch when its matching external scale is exactly zero. When the external scale is nonzero, only the retained relative gate applies.

Let

[ \Delta\mathbf q_{inc} =\mathbf q_{trial}-\mathbf q_{step-start}. ]

Correction channels use

[ \eta_\Delta^u= \frac{|\Delta\mathbf u_f|\infty} {|\Delta\mathbf u{inc,f}|\infty}, \qquad \eta\Delta^\phi= \frac{|\Delta\boldsymbol\phi_f|\infty} {|\Delta\boldsymbol\phi{inc,f}|_\infty} ]

for nonzero denominators and pass iff each is (\le10^{-2}). An exact-zero denominator instead passes iff the fresh numerator is at most (10^{-3}) length or (10^{-3}) dimensionless rotation. Empty free channels have exact-zero numerator/denominator and pass their typed zero-scale rule.

All four gates are combined by logical AND. Constrained coordinates are excluded. These are approved solver-convergence values, not external-reference tolerances.

15. Verification Identities

The following are formulation identities whose numerical norms/tolerances are set by Numerical Review.

15.1 Linear limit

At (\mathbf q=\mathbf0), stress-free reference:

[ \mathbf A_0=\mathbf T_p,\qquad \mathbf r_{20}^{phys}=\mathbf0,\qquad \mathbf K_{20}^{geo}=\mathbf0, ]

[ \mathbf K_{phys}^{24} =\mathbf T_p^T\mathbf K_{20}^{linear}\mathbf T_p. ]

The nonlinear Green–Lagrange/MITC physical strain, stress, resultants and recovery reduce to the approved linear MITC4 signs, component orders and 2 x 2 x 2 rule. A pure drill basis reduces to (k_d\mathbf I_4). Objective corotation adds only the coupling required to make finite rigid motion drill-free and contributes neither physical recovery nor stress.

15.2 Objectivity and six rigid modes

For any (\mathbf R\in SO(3)), (\mathbf c\in\mathbb R^3),

[ \mathbf x^+=\mathbf R\mathbf x+\mathbf c,\qquad \mathbf Q_I^+=\mathbf R\mathbf Q_I ]

implies

[ \overline{\mathbf E}^+=\overline{\mathbf E},\quad \Pi_{phys}^+=\Pi_{phys},\quad \boldsymbol\gamma_e^+=\boldsymbol\gamma_e,\quad \Pi_{drill}^+=\Pi_{drill}. ]

At a stress-free reference with no loads, three translations and three proper rigid rotations produce zero residual and zero energy. After objective drill regularization, an isolated valid element has exactly these six physical rigid modes; pure numerical drill is not an additional zero mode.

15.3 Tangent consistency

For a finite valid trial and deterministic direction (\mathbf p),

[ \boxed{ \mathbf K_T(\mathbf q)\mathbf p =D\mathbf r(\mathbf q)[\mathbf p].} ]

The directional comparison must exercise separately and together:

• material term;
• Green–Lagrange geometric term;
• (\mathbf A^T\mathbf K_{20}\mathbf A);
• (\sum r_{20,a}\mathbf H_a^\Phi);
• fixed-spatial source-moment projector, rotation-coordinate and map-Hessian terms;
• objective drill gradient and exact Hessian.

The two perturbed states share the same fixed reference, load factor and committed base. Perturbation schedule and tolerance are Numerical Review ownership.

15.4 Symmetry, equilibrium and state

• With no applied nodal moment, the complete unreduced physical/drill tangent is symmetric in exact arithmetic. With a fixed-spatial nodal moment, the exact dead-moment work tangent is generally nonsymmetric and must not be symmetrized.
• Element and assembled internal actions satisfy virtual-work covariance under stable scatter/assembly.
• Accepted free residual is equilibrium evidence; constrained full residual is reaction.
• Global resultant force and moment identities in §13.4 include nodal couples and current force arms.
• Every committed factor has one accepted frame; no failed trial changes the preceding (\mathbf u,\boldsymbol\phi,\mathbf Q,\mathbf d), force, reaction or recovery identity.

15.5 Patch and recovery checks

• partition of unity, Kronecker delta and derivative-sum identities;
• constant membrane, bending, shear and twist continuation from the linear limit;
• direct/tied shear component identity and sign;
• bottom/middle/top stress and section-resultant thickness-moment checks;
• regular, rigidly rotated and admissibly distorted patches;
• no physical strain/stress/resultant change from pure numerical drill;
• stable source, element, Gauss, section and component ordering.

16. Numerical Risks

risk	consequence	formulation control / downstream evidence
Rigid-body mode loss	artificial force/energy	objective physical and drill potentials; exactly six physical rigid modes
Patch-test failure	inconsistent interpolation/sign	direct membrane/bending plus identical MITC projection of strain and two variations
Tangent asymmetry/inconsistency	degraded or failed Full Newton	exact material, geometric, map, drill Hessians and exact dead-moment work derivative; directional derivative review
Indefinite tangent	valid near-buckling state or numerical failure	do not require positive definiteness globally; distinguish expected physical indefiniteness from singular/malformed tangent
Hourglass / zero-energy mode	singular element	no reduced integration; 2 x 2 x 2 full rule; objective drill removes only numerical drill modes
Transverse-shear locking	overly stiff thin response	canonical MITC tying; thickness/mesh characterization
Membrane locking	overly stiff curved/distorted original MITC4	documented known limitation; do not claim MITC4+ performance
Volumetric locking	not applicable to approved plane-stress shell	no 3-D incompressible claim; material range remains (\nu<0.5)
Director/swing antipode	non-unique physical chart	fail closed at (d^0\cdot d=-1); no silent axis choice
Exponential/log branch	non-unique global UR target/state	principal (|\phi|<\pi); downstream branch tolerance
Drill twist branch	discontinuous numerical potential	fail closed at non-unique (\pm\pi) twist
Nonobjective drilling	artificial rigid energy	polar element corotation and exact potential Hessian
State-dependent drill scale	missing scale derivative/nondeterminism	freeze initial approved (k_d); no trial recomputation
Director-parallel applied moment	physical work mixed with drill	immutable-source current-trial fail-closed ratio; exact physical-work projection; no drill channel
Moment derivative omission	inconsistent residual/tangent	immutable (M_{src}^\star), exact projector/coordinate/map derivative, no follower direction
Reversed/singular Jacobian	invalid basis, sign or inverted shell	(J_0>0), (j>0), (\det F_c>0) at every required location
Distortion/warpage	conditioning and membrane locking	exact point inventory; no omitted/averaged failed location
Nonfinite trial	corrupted state/output	reject before solve acceptance or commit
Constraint/target confusion	wrong load path or reaction	absolute path applied before iteration; constrained Newton correction exactly zero
Mixed units in convergence	false pass/fail	force/moment and translation/rotation channels remain separate
Determinism loss	irreproducible reference/evidence	stable source order, local buffers, fixed reduction and diagnostic order

17. Requirement and Research Traceability

requirement / gap	formulation closure
FESA-REQ-NSNR-003	§§2, 6–10: original-style MITC4, homogeneous isotropic materially linear TL shell
FESA-REQ-NSNR-041	§§2, 6–9, 12: one fixed initial reference, Green–Lagrange / second-Piola pair
FESA-REQ-NSNR-074	§3: 6 global DOF, dimensionless correction, committed/trial (Q,d), multiplicative update
FESA-REQ-NSNR-006	§§11, 14: complete same-trial residual/tangent and fresh Full Newton solve
FESA-REQ-NSNR-007/069/070	§14.1: four typed inclusive gates and constrained exclusion
FESA-REQ-NSNR-008/014	§14: finite data, iterations 1..100, unchanged correction, no fallback
FESA-REQ-NSNR-046/048	§§4, 14: one proportional factor and fixed-step lifecycle
FESA-REQ-NSNR-018	§§4.4, 5, 11.2: immutable fixed-global source coefficients, exact tangent-space work covector/derivative and no follower direction
FESA-REQ-NSNR-063	§§4.2–4.3: factor-one absolute target and zero constrained Newton correction
FESA-REQ-NSNR-011	§§3.5, 14: separate trial and atomic complete commit/rollback
FESA-REQ-NSNR-019	§§11.4, 13.1, 13.4, 14: deterministic full residual, dimensioned reaction covectors and force/moment balance
FESA-REQ-NSNR-030	§15: linear limit, tangent, equilibrium, objectivity and rollback identities
O-ROT-01	§§3.2–3.4: fixed-reference exponential coordinate; exact left-spatial multiplicative update and derivatives
O-ROT-02	§§3, 6–8, 12: nodal triad/director, exact director interpolation, point frame/branch validation
O-MAP-01	§§7, 11: exact (\Phi:\mathbb R^{24}\to\mathbb R^{20}), (A), (H^\Phi), complete map curvature
O-DRILL-01	§10: objective polar-corotated twist potential, frozen scale, exact residual/Hessian, six rigid modes
O-DRILL-02	§§4.4, 5, 10, 11.2: current-director source admissibility, exact physical-work projection and no drill load channel
O-MOM-01	§§5, 11.2–11.3: immutable fixed-spatial source, projected physical swing covector, exact work derivative and no follower direction
O-BC-01	§§3.2–3.3, 4.2–4.3, 14: absolute component path, exact multiplicative triad update, (\Delta q_c=0)

Theory basis follows research sources S1/S6 for fixed-reference nonlinear continuum-based MITC shell and tying/integration, S2–S5 for finite rotation manifold and exact potential derivatives, and S7 for same-state Full Newton. Local wiki cross-checks were Total Lagrangian Shell Formulation, Green-Lagrange Strain Linearization, MITC4 Shell Element, MITC Shell Kinematics, Continuum Mechanics Based Four-Node Shell Element, Assumed Transverse Shear Strain Interpolation, and Nonlinear Finite Element Analysis. Project requirements and approved repository contracts take precedence over those navigation pages.

18. Open Issues and Downstream Handoff

18.1 Formulation closures

O-ROT-01, O-ROT-02, O-MAP-01, O-DRILL-01, O-DRILL-02, O-MOM-01 and O-BC-01 are closed by the explicit equations and fail-closed boundaries in this document. No placeholder derivation remains for implementation of their residual or tangent.

18.2 Numerical Review-owned decisions

The following remain intentionally open to Numerical Review and do not alter this formulation:

• O-TOL-01: scale-aware tolerances for tangent derivative, symmetry, objectivity, triad invariants, equilibrium, branch proximity and Jacobian validity;
• O-REF-01/O-REF-02: exact stored nonlinear input/CSV and independent-source readiness;
• O-BENCH-01: exact blocking large-response benchmark and quantities;
• O-STEP-01: endpoint arithmetic and fixed increment adequacy without cutback;
• O-MMS-01: whether shell MMS remains nonblocking.

Numerical Review must independently check dimensions, signs, exact derivative identities, the polar-corotation drill construction, six rigid modes, prescribed rotation path, moment rejection, tangent symmetry/consistency, locking limitations, and benchmark traversability. This status is ready-for-numerical-review, not numerically approved and not a downstream implementation authorization.

18.3 Boundary of later stages

• Exact HDF5 datasets, log representation, parser/semantic validation and public CLI behavior remain I/O Definition ownership.
• Exact reference artifacts, comparison families and tolerance application remain Numerical Review/reference contract ownership under docs/TOLERANCE.md.
• C++ type/API/file ownership, sparse implementation and tests remain later planning and implementation ownership.
• No Abaqus, Nastran or other reference solver execution, artifact creation, artifact modification, or internal-equivalence claim is authorized by this document.