FESADev/docs/MITC4_FORMULATION.md

MITC4 Formulation

Purpose

This document defines the baseline MITC4 formulation target for FESA Phase 1.

It is intentionally a formulation contract, not implementation code. Exact formulas should be added and reviewed before coding the MITC4 element.

Source Basis

• Dvorkin and Bathe's four-node shell element paper presents a continuum-mechanics-based, non-flat, general quadrilateral shell element for thin and thick shells and nonlinear analysis: https://web.mit.edu/kjb/www/Publications_Prior_to_1998/A_Continuum_Mechanics_Based_Four-Node_Shell_Element_for_General_Nonlinear_Analysis.pdf
• The paper identifies transverse shear locking as a key problem in simple 4-node shell interpolation and motivates modified transverse shear treatment: https://web.mit.edu/kjb/www/Publications_Prior_to_1998/A_Continuum_Mechanics_Based_Four-Node_Shell_Element_for_General_Nonlinear_Analysis.pdf
• OpenSees describes ShellMITC4 as a bilinear isoparametric shell element with modified shear interpolation, four counter-clockwise nodes, and six DOFs per node: https://opensees.berkeley.edu/wiki/index.php/Shell_Element
• The MITC benchmark paper states that the MITC method is used to remedy shell locking and that the standard MITC4 employs MITC treatment for transverse shear strains; it also notes that Abaqus S4 uses Dvorkin-Bathe transverse shear interpolation: https://web.mit.edu/kjb/www/Principal_Publications/Performance_of_the_MITC3%2B_and_MITC4%2B_shell_elements_in_widely_used_benchmark_problems.pdf
• Abaqus finite-strain shell theory documentation provides useful comparison context for S4/S4R geometry, interpolation, orientation update, and transverse shear treatment, but FESA Phase 1 is linear static: https://abaqus-docs.mit.edu/2017/English/SIMACAETHERefMap/simathe-c-finitestrainshells.htm

Phase 1 Target

Phase 1 implements a clear MITC4 baseline formulation and passes reference benchmarks before performance optimization.

Scope:

• 4-node quadrilateral shell.
• Linear static analysis.
• Linear isotropic elastic material.
• Homogeneous shell section.
• 6 DOFs per node.
• Small-strain formulation for Phase 1.
• Transverse shear interpolation based on MITC4 assumptions.
• Abaqus-compatible result signs.

Non-scope:

• S4R reduced-integration behavior.
• Hourglass control.
• Composite sections.
• Material nonlinearity.
• Geometric nonlinearity.
• Pressure loads.
• Thermal-stress coupling.
• Mesh quality diagnostics.

Nodal DOFs

Each node has:

UX, UY, UZ, RX, RY, RZ

Rules:

• Translational DOFs are global translations.
• Rotational DOFs are rotations about global or transformed axes following Abaqus component convention.
• RZ is retained as a drilling DOF.
• Drilling stiffness is artificial in Phase 1 and must be parameterized.

Element Input Contract

MITC4Element requires:

• four node ids in Abaqus S4 order.
• four node coordinates.
• shell thickness.
• linear elastic material constants E and nu.
• drilling stiffness parameter.
• element id and property id for diagnostics and output.

Node ordering:

• Input node order follows Abaqus S4 convention.
• Positive normal follows the right-hand rule around the nodes.
• FESA maps Abaqus TYPE=S4 to MITC4.
• Abaqus TYPE=S4R is not supported in Phase 1.

Coordinate Frames

The exact local basis construction must be completed before MITC4 implementation.

Minimum requirements:

• Define a local shell normal from the quadrilateral geometry.
• Define local in-plane axes e1 and e2 so that e1, e2, and normal form a right-handed basis.
• Preserve Abaqus-compatible output signs.
• Document behavior for non-planar quadrilaterals.
• Use the same convention consistently for stiffness, stress/strain recovery, and result output.

Recommended Phase 1 convention:

• Use the element midsurface geometry to compute an average normal.
• Use a projected global axis to define the local 1-direction when possible, matching Abaqus convention conceptually.
• Fall back to a stable element-edge-based direction when the projected global axis is nearly parallel to the normal.
• Record the final algorithm in this document before coding.

Shape Functions

Baseline quadrilateral bilinear interpolation:

N1 = 0.25 * (1 - r) * (1 - s)
N2 = 0.25 * (1 + r) * (1 - s)
N3 = 0.25 * (1 + r) * (1 + s)
N4 = 0.25 * (1 - r) * (1 + s)

where r, s are natural coordinates in [-1, 1].

Implementation requirements:

• Compute shape function derivatives with respect to natural coordinates.
• Build the surface Jacobian and local derivatives.
• Detect invalid or near-zero Jacobian as a singular/invalid element diagnostic, not as a mesh quality metric.

Strain Treatment

The baseline element separates:

• membrane strain terms.
• bending curvature terms.
• transverse shear strain terms.
• artificial drilling stabilization.

MITC4 requirement:

• Use standard displacement interpolation for membrane and bending terms in Phase 1.
• Use MITC transverse shear interpolation to alleviate shear locking.
• Do not replace MITC4 with plain full-integration Reissner-Mindlin Q4.

The exact tying point equations and shear interpolation formula must be added from the selected primary source before implementation.

Numerical Integration

Initial Phase 1 plan:

• In-plane integration: 2x2 Gauss for membrane/bending/shear stiffness unless the final formulation requires a different split.
• Thickness integration: homogeneous linear elastic section may be integrated analytically or with a documented simple rule.
• Benchmark literature commonly reports 2x2 in-plane Gauss integration for S4/MITC4-style 4-node elements and 2-point thickness integration in comparative shell studies.

Rules:

• Do not introduce reduced integration or hourglass control for S4R behavior in Phase 1.
• Do not optimize integration before reference benchmarks pass.
• Integration point ordering for output must be documented before stress/strain reference comparisons.

Drilling DOF Stabilization

Decision:

• Phase 1 uses small artificial drilling stiffness.

Requirements:

• Expose a parameter such as drilling_stiffness_scale.
• Provide a deterministic default.
• Make the default small enough not to dominate physical response.
• Include benchmark sensitivity checks if reference results are sensitive to the value.
• Report the value in result metadata.

Open default proposal:

k_drill = alpha * representative_element_stiffness

where alpha should be selected only after reviewing formulation sources and early reference cases.

Element Outputs

Phase 1 minimum:

• element stiffness matrix.
• element equivalent nodal internal force for full-vector residual/reaction recovery.
• optional stress/strain output after displacement benchmarks are stable.

Future output:

• local shell stresses.
• local shell strains.
• section forces and moments.
• integration point and section point data.

Required Element-Level Tests

Before integration with the global solver:

• shape functions sum to one.
• derivatives satisfy expected bilinear identities.
• element stiffness dimensions are 24 x 24.
• stiffness is symmetric for linear elastic Phase 1.
• rigid body translations produce near-zero internal strain energy.
• rigid body rotations do not create physical membrane/bending stiffness beyond documented drilling effects.
• constant membrane patch behavior.
• bending-dominated sanity case.
• drilling stiffness sensitivity check.

Reference Benchmarks

MITC4 baseline acceptance should include:

• single-element membrane test.
• single-element bending test.
• cantilever shell strip.
• simply supported square plate.
• Scordelis-Lo roof.
• pinched cylinder.
• twisted beam.

Distorted mesh tests should be added after the baseline passes, but Phase 1 does not implement general mesh quality diagnostics.

Future Extensions

Geometric nonlinearity:

• Add updated geometry, current frame handling, tangent stiffness, and Newton-Raphson integration.
• Preserve AnalysisState element/internal state extension points.

Thermal-stress coupling:

• Add temperature field state.
• Add thermal strain contribution.
• Add material expansion data.
• Add result fields for temperature and thermal strain/stress.

S4R:

• Add only after a separate ADR and formulation document update.
• Requires reduced integration and hourglass control decisions.

Open Decisions Before Coding

• Exact MITC4 transverse shear tying point formula.
• Exact element local basis for warped quads.
• Exact drilling stiffness default.
• Exact stress/strain recovery locations.
• Whether Phase 1 reports S, E, and SF, or only U and RF.
• Whether local coordinate transforms from Abaqus input are deferred or rejected.