FESADev/docs/MITC4_FORMULATION.md

MITC4 Formulation

Purpose

This document is the formulation contract for FESA's Phase 1 four-node MITC4 shell element.

It supersedes the earlier simplified Phase 1 baseline that used an averaged-edge shell basis, analytic thickness integration, and drilling_stiffness_scale = 1.0e-6. The next Phase 1 implementation pass must treat this document as the MITC4 gate before coding or evaluating shell stiffness.

The goal is not to reproduce Abaqus internals. The goal is to implement a documented, testable Dvorkin-Bathe-style MITC4 element that can be compared against stored Abaqus S4 reference artifacts without running Abaqus locally.

Local Source Basis

The formulation below is based on the local papers in docs/Paper/:

• docs/Paper/AContinuumMechanicsBasedFourNodeShellElementforGeneralNonlinearAnalysis/
  • Dvorkin and Bathe's original four-node non-flat shell element.
  • Defines the degenerated-continuum geometry, displacement interpolation, convected transverse shear interpolation, constitutive transformation, and benchmark expectations.
• docs/Paper/FourNodeQuadrilateralShellElementMITC4/
  • Restates the MITC4 geometry and displacement interpolation, gives explicit midside shear tying relations, and reports patch-test and Scordelis-Lo verification behavior.
• docs/Paper/유한요소해석법을이용한쉘구조물의동적좌굴해석/
  • Provides the Korean thesis derivation for covariant strain, plane-stress material matrix with shear correction, 6-DOF shell transformation, drilling stabilization, and nonlinear extension notes.
• docs/Paper/2007쉘구조물의유한요소해석에대하여/
  • Provides shell-model background, locking behavior, asymptotic behavior, and benchmark/test interpretation.
• docs/Paper/mitc공부/
  • Treated as supporting study notes only. Use the papers above for binding formulas.

Phase 1 Scope

Phase 1 implements:

• Abaqus TYPE=S4 mapped to FESA MITC4.
• Four-node quadrilateral shell in Abaqus S4 node order.
• Linear static, small-strain, linear isotropic elastic analysis.
• Homogeneous shell section with one thickness value per element for the initial implementation.
• Six global DOFs per node: UX, UY, UZ, RX, RY, RZ.
• MITC transverse shear interpolation in a convected coordinate system.
• Constrained DOF elimination, full-vector reaction recovery, and minimum output fields U and RF.

Phase 1 does not implement:

• Abaqus S4R, reduced integration, or hourglass control.
• Geometric nonlinearity, material nonlinearity, dynamics, pressure loads, thermal-stress coupling, composite sections, or mesh quality diagnostics.
• Abaqus *Orientation, user-defined nodal normals, or shell offset behavior unless added by a later ADR.

Natural Coordinates and Node Order

FESA uses the standard bilinear quadrilateral natural coordinates:

node 1: (xi, eta) = (-1, -1)
node 2: (xi, eta) = (+1, -1)
node 3: (xi, eta) = (+1, +1)
node 4: (xi, eta) = (-1, +1)
zeta in [-1, +1] through the thickness

Shape functions:

N1 = 0.25 * (1 - xi) * (1 - eta)
N2 = 0.25 * (1 + xi) * (1 - eta)
N3 = 0.25 * (1 + xi) * (1 + eta)
N4 = 0.25 * (1 - xi) * (1 + eta)

Midside MITC tying points in the FESA convention:

A = ( 0, -1)  edge 1-2
B = (-1,  0)  edge 1-4
C = ( 0, +1)  edge 4-3
D = (+1,  0)  edge 2-3

Dvorkin-Bathe papers may use a natural-coordinate orientation where the signs on the A/C interpolation are reversed relative to the FESA convention above. FESA implementation must preserve the edge definitions above and use the FESA interpolation formulas in the MITC shear section.

Degenerated-Continuum Geometry

The shell is represented as a 3D continuum degenerated through the thickness:

X(xi, eta, zeta) =
  sum_k N_k(xi, eta) X_k
  + zeta / 2 * sum_k t_k N_k(xi, eta) Vn_k

where:

• X_k is the midsurface coordinate of node k.
• t_k is the nodal thickness. Phase 1 uses t_k = t for all four nodes.
• Vn_k is the nodal director. The papers allow the director to be non-normal to the midsurface; Phase 1 derives it from geometry because Abaqus nodal normals are not in the parser subset.

Phase 1 director policy:

• If no parsed nodal directors exist, use the element-center midsurface normal for all four nodal directors:

G1_c = dX_mid / dxi  at (0, 0)
G2_c = dX_mid / deta at (0, 0)
Vn_k = normalize(G1_c x G2_c), k = 1..4

• If G1_c x G2_c is near zero, reject the element with an invalid/singular element diagnostic.
• This flat-or-mildly-warped policy is sufficient for the Phase 1 S4 baseline. A true non-flat nodal-director policy requires a later ADR and benchmark update.

The nodal local director axes are:

V1_k = normalize(EY x Vn_k)
V2_k = Vn_k x V1_k

If EY x Vn_k is near zero, use a deterministic fallback axis before normalization:

V1_k = normalize(EZ x Vn_k)
if still near zero, V1_k = normalize(EX x Vn_k)
V2_k = Vn_k x V1_k

The implementation must keep V1_k, V2_k, and Vn_k orthonormal and right-handed.

Displacement and Rotation Interpolation

The local five-DOF MITC4 displacement field is:

u(xi, eta, zeta) =
  sum_k N_k u_k
  + zeta / 2 * sum_k t_k N_k q_k

q_k = -V2_k * alpha_k + V1_k * beta_k

where:

• u_k is the nodal translation vector.
• alpha_k is the rotation about V1_k.
• beta_k is the rotation about V2_k.
• Drilling rotation gamma_k about Vn_k does not enter the continuum strain field.

FESA stores six global DOFs per node. Global nodal rotations are mapped to local director rotations by:

theta_k = [RX_k, RY_k, RZ_k]^T

[alpha_k, beta_k, gamma_k]^T = L_k theta_k

L_k =
  [ EX dot V1_k   EY dot V1_k   EZ dot V1_k
    EX dot V2_k   EY dot V2_k   EZ dot V2_k
    EX dot Vn_k   EY dot Vn_k   EZ dot Vn_k ]

The element-level local-to-global transformation is block diagonal:

u_local_e = T_e u_global_e
K_global_e = T_e^T K_local_e T_e
f_global_e = T_e^T f_local_e

Each node block maps [UX, UY, UZ, RX, RY, RZ] to [u_x, u_y, u_z, alpha, beta, gamma].

Covariant Strain Definition

Use convected covariant components for the MITC4 strain construction:

g_i = dX / dr_i, where r_i = xi, eta, zeta

eps_ij =
  0.5 * (du/dr_i dot g_j + du/dr_j dot g_i)

For the Phase 1 small-strain implementation:

• Membrane and bending terms are derived directly from the degenerated-continuum displacement field.
• eps_33 = 0 is enforced by the shell assumption.
• The transverse shear rows eps_13 and eps_23 are not taken directly at the Gauss point. They are replaced by the MITC assumed transverse shear interpolation below.

Use the internal strain vector order:

[ eps_11, eps_22, eps_33, gamma_23, gamma_13, gamma_12 ]^T

where engineering shear strains are:

gamma_ij = 2 * eps_ij

If an implementation uses another row order internally, it must include an explicit permutation test against this documented order.

MITC Transverse Shear Interpolation

The key MITC4 assumption is that the transverse shear tensor components are interpolated from midside tying points, instead of being evaluated directly at every Gauss point from the displacement field.

The original paper form is:

eps_13 = 0.5 * (1 + r2) * eps_13^A + 0.5 * (1 - r2) * eps_13^C
eps_23 = 0.5 * (1 + r1) * eps_23^D + 0.5 * (1 - r1) * eps_23^B

With FESA's Abaqus-style natural coordinate convention and tying-point labels, use:

eps_13_MITC(xi, eta) =
  0.5 * (1 - eta) * eps_13^A
  + 0.5 * (1 + eta) * eps_13^C

eps_23_MITC(xi, eta) =
  0.5 * (1 - xi) * eps_23^B
  + 0.5 * (1 + xi) * eps_23^D

The tying values are the direct covariant shear strains evaluated at midside points:

eps_13^A = eps_13_DIRECT( 0, -1)
eps_13^C = eps_13_DIRECT( 0, +1)
eps_23^B = eps_23_DIRECT(-1,  0)
eps_23^D = eps_23_DIRECT(+1,  0)

Implementation rule:

• Build the direct covariant strain-displacement rows from the geometry and displacement interpolation.
• Evaluate only the needed direct transverse shear rows at A, B, C, and D.
• Interpolate those rows to the integration point using the formulas above.
• Do not compute MITC shear from local Cartesian gamma_xz/gamma_yz first; the tying is on convected tensor components.

For sign and regression checks, define:

q_k = -V2_k * alpha_k + V1_k * beta_k

The paper edge formulas can be used as a diagnostic check for the direct tying rows:

eps_13^A ~ 1/8 * [
  (u1 - u2) dot 0.5 * (t1 Vn1 + t2 Vn2)
  + (X1 - X2) dot 0.5 * (t1 q1 + t2 q2)
]

eps_13^C ~ 1/8 * [
  (u4 - u3) dot 0.5 * (t4 Vn4 + t3 Vn3)
  + (X4 - X3) dot 0.5 * (t4 q4 + t3 q3)
]

eps_23^B ~ 1/8 * [
  (u1 - u4) dot 0.5 * (t1 Vn1 + t4 Vn4)
  + (X1 - X4) dot 0.5 * (t1 q1 + t4 q4)
]

eps_23^D ~ 1/8 * [
  (u2 - u3) dot 0.5 * (t2 Vn2 + t3 Vn3)
  + (X2 - X3) dot 0.5 * (t2 q2 + t3 q3)
]

These compact formulas are not a license to skip the general direct-row evaluation. They are a sign/orientation cross-check derived from the same midside strain definitions.

Local Cartesian Frame and Material Law

The constitutive law is defined in a local orthonormal Cartesian frame at each integration point with sigma_33 = 0.

Build the integration-point local basis from the covariant basis:

e3_hat = normalize(g_3)
e1_hat = normalize(g_2 x e3_hat)
e2_hat = e3_hat x e1_hat

If g_2 x e3_hat is near zero, use a deterministic fallback that still creates a right-handed orthonormal basis and emits a diagnostic if no valid basis exists.

For isotropic linear elasticity, use the local plane-stress shell material matrix with transverse shear correction kappa = 5/6:

D_hat = E / (1 - nu^2) *
  [ 1   nu  0   0                    0                    0
    nu  1   0   0                    0                    0
    0   0   0   0                    0                    0
    0   0   0   kappa*(1-nu)/2       0                    0
    0   0   0   0                    kappa*(1-nu)/2       0
    0   0   0   0                    0                    (1-nu)/2 ]

This matrix uses the documented strain order:

[ eps_11, eps_22, eps_33, gamma_23, gamma_13, gamma_12 ]^T

For a general non-flat element, transform the local Cartesian material relation to the convected coordinate strain vector:

E_hat = T_xi_to_x E_convected
D_convected = T_xi_to_x^T D_hat T_xi_to_x

The transformation matrix must be constructed from direction cosines between the local orthonormal basis and contravariant/covariant convected bases. For Phase 1, a flat-element implementation may compute the B-matrix directly in the local Cartesian frame only if tests prove it is equivalent to the convected formulation for the accepted reference cases.

Element Stiffness and Numerical Integration

The baseline stiffness is:

K_e = integral_V B^T D B dV

Phase 1 baseline integration:

• Use 2 x 2 x 2 Gauss integration in (xi, eta, zeta) for the first rebuilt MITC4 implementation.
• This follows the thesis derivation's 2 x 2 x 2 natural-coordinate integration and the Dvorkin-Bathe elastic recommendation of 2 x 2 in the midsurface with two through-thickness points.
• Do not introduce reduced integration or hourglass control.
• Analytic through-thickness integration may be introduced later only after equivalence tests preserve the MITC shear, drilling, and reference displacement behavior.

At every integration point:

1. Compute g_i, the metric/Jacobian, local Cartesian basis, and material transform.
2. Build direct membrane/bending strain rows.
3. Replace transverse shear rows with MITC interpolated rows.
4. Accumulate B^T D B |J| w_xi w_eta w_zeta.

Detect near-zero Jacobian, invalid basis, missing property/material, or non-positive thickness as invalid/singular element diagnostics. Do not report these as mesh quality diagnostics in Phase 1.

Drilling DOF Stabilization

The continuum MITC4 formulation has five physical DOFs per node. FESA keeps six global DOFs, so the local drilling rotation gamma has no physical strain contribution and must be weakly stabilized.

Baseline decision:

drilling_stiffness_scale = 1.0e-3
k_ref = min_positive_diagonal(K_local_without_drilling)
k_drill = drilling_stiffness_scale * k_ref

Rules:

• K_local_without_drilling is the element local stiffness after MITC integration and before artificial drilling stabilization.
• Use the minimum positive diagonal entry over the physical local DOFs. Ignore zero, negative, NaN, or non-finite entries.
• Add k_drill to each nodal local gamma diagonal before transforming to global coordinates.
• If no positive reference diagonal exists, emit a singular/invalid element diagnostic.
• The scale must be configurable and written to result metadata or analysis diagnostics.
• Reference comparisons must include a drilling sensitivity check if displacement results change materially with the scale.

This replaces the earlier arbitrary 1.0e-6 * E * thickness rule. The 1.0e-3 scale follows the local thesis derivation's six-DOF stabilization note and is still an artificial numerical parameter, not a physical stiffness.

Internal Force and Reaction Recovery

For Phase 1 linear static analysis:

f_int_e = K_e u_e
R_full = K_full U_full - F_full

Rules:

• Reactions must be recovered from the full global vector, not only from the reduced free-DOF system.
• RF output follows the Abaqus-style component convention documented in docs/NUMERICAL_CONVENTIONS.md.
• Element stress/strain/resultant outputs are not mandatory until displacement and reaction reference behavior is stable.

Required Tests Before Reimplementation

The rebuilt MITC4 element must be protected by tests before implementation code is accepted.

Element-level tests:

• Shape functions sum to one and have correct derivatives at center, corners, and tying points.
• Center normal, local director axes, and integration local Cartesian axes are orthonormal and right-handed.
• The local-to-global rotation transform maps global rotations to alpha, beta, gamma as documented.
• Direct tying rows at A, B, C, D agree with finite-difference strain perturbations.
• MITC shear rows at Gauss points are interpolated from tying rows, not evaluated directly at the Gauss points.
• Stiffness is symmetric for linear elastic Phase 1.
• Rigid-body translations and rotations produce near-zero physical strain energy, with only documented drilling stabilization effects.
• The element has no spurious zero-energy modes beyond expected rigid-body behavior after drilling stabilization is accounted for.

Patch and benchmark tests:

• Constant membrane states.
• Pure bending.
• Pure shear.
• Pure twist.
• Thin cantilever or plate strip to expose shear locking.
• Scordelis-Lo roof after the basic flat tests pass.
• Pinched cylinder or twisted beam as later Phase 1/Phase 2 benchmark expansion.

Reference regression tests:

• Compare FESA U against stored Abaqus displacement CSV files using absolute and relative tolerances.
• Keep RF verified by full-vector equilibrium until matching Abaqus reaction CSV files are provided or explicitly deferred.

Current Reference Artifacts

Stored artifacts currently known:

• references/quad_01.inp
• references/quad_01_displacements.csv
• references/quad_02.inp
• references/quad_02_displacements.csv
• references/quad_02_reactionforces.csv

Compatibility notes:

• quad_01.inp contains S4R, Part/Assembly/Instance, *Density, and NLGEOM=YES; it remains future compatibility provenance, not a Phase 1 parser acceptance case.
• quad_02.inp contains TYPE=S4, which is the correct element target for Phase 1, but it also uses Part/Assembly/Instance. The normalized quad_02_phase1.inp remains the executable Phase 1 input path. Do not silently expand parser support while implementing or verifying the element.
• quad_02_reactionforces.csv is now available for Abaqus RF/RM comparison. The current node-wise RF comparison does not pass; keep this as a formulation/solver verification gap until the mismatch is explained or fixed.

Implementation Checklist

The next MITC4 implementation pass should proceed in this order:

1. Lock the shape-function, natural-coordinate, node-order, and tying-point tests.
2. Implement geometry/director/local-axis construction with diagnostics.
3. Implement the five-DOF degenerated-continuum displacement interpolation.
4. Implement global six-DOF to local [alpha, beta, gamma] transformation.
5. Implement direct covariant strain rows and finite-difference tests.
6. Replace transverse shear rows with MITC tying interpolation.
7. Implement material transform and 2 x 2 x 2 stiffness integration.
8. Add drilling stabilization in local coordinates.
9. Transform element stiffness and internal force to global coordinates.
10. Run element patch tests before assembly/reference regression work.

Future Extensions

Geometric nonlinearity:

• Update current geometry and nodal directors in AnalysisState.
• Add finite-rotation updates for V1, V2, and Vn.
• Add tangent stiffness, geometric stiffness, and Newton-Raphson convergence checks.

Thermal-stress coupling:

• Add temperature field state and thermal strain contribution.
• Add material expansion data and thermal output fields.

S4R:

• Add only after a separate ADR and formulation update.
• Requires reduced integration, hourglass control, and separate reference artifacts.

Stress/resultant output:

• Define section point locations, component labels, local coordinate convention, and Abaqus CSV comparison format before making S, E, or SF mandatory.

Remaining Open Decisions

The MITC4 stiffness formulation is now defined for the Phase 1 rebuild. Remaining decisions are outside the first stiffness pass:

• Exact stress/strain/resultant recovery locations and component labels.
• Whether to parse Abaqus *Orientation or nodal shell normals.
• Whether quad_02.inp should be normalized into the existing Phase 1 parser subset or used to justify a dedicated Part/Assembly/Instance parser sprint.
• Reference tolerances for quad_02 after the input compatibility path is chosen.