AbaqusSubroutineDev/docs/formulations/uel-3d-euler-beam.md

3D Euler-Bernoulli Beam UEL Formulation

Metadata

• feature_id: uel-3d-euler-beam
• source_requirement: docs/requirements/uel-3d-euler-beam.md
• source_research: docs/research/uel-3d-euler-beam-research.md
• status: ready-for-numerical-review
• owner_agent: formulation-agent
• date: 2026-06-11

Scope and Assumptions

• analysis_type: linear static small-displacement
• Abaqus entry point: UEL
• element_type: two-node 3D Euler-Bernoulli beam/frame element
• topology: node 1 to node 2 defines the positive local 1 axis
• deformation: small displacement and small rotation only
• material_model_boundary: linear elastic section stiffness with axial, Saint-Venant torsion, and two principal bending stiffnesses
• coordinate_system: Abaqus global nodal coordinates for AMATRX and RHS; element local coordinates only for intermediate beam calculations
• units: user-consistent units; no unit conversion inside the element
• loads: no distributed, body, thermal, follower, or element-generated loads in first scope
• dynamics: no mass, damping, inertia, modal, or transient contribution in first scope

This document defines the mathematical element contract. It does not define the final Abaqus ABI signature, PROPS ordering, JPROPS ordering, source file layout, or Fortran implementation plan.

Primary Variables and DOFs

The element has two nodes and 12 active DOFs. Global and local element vectors use the same component ordering; the difference is only the coordinate basis.

q = [
  u1, v1, w1, th1_1, th2_1, th3_1,
  u2, v2, w2, th1_2, th2_2, th3_2
]^T

Where:

• u is translation along local 1.
• v is translation along local 2.
• w is translation along local 3.
• th1 is rotation about local 1.
• th2 is rotation about local 2.
• th3 is rotation about local 3.

The Abaqus-facing DOF order remains node 1 U1, U2, U3, UR1, UR2, UR3, followed by node 2 U1, U2, U3, UR1, UR2, UR3. In the local vector above, each translational or rotational component is the projection of the corresponding Abaqus global vector onto the local frame.

Euler-Bernoulli bending imposes the section-normal constraints:

th3 =  dv/dx
th2 = -dw/dx

The sign difference follows the right-handed local 1-2-3 basis: a positive rotation about local 3 changes the centerline tangent in the local 2 direction, while a positive rotation about local 2 changes it in the negative local 3 direction.

Required Material and Section Parameters

The first-scope formulation requires finite, strictly positive:

parameter	meaning	stiffness role
E	Young's modulus	bending and axial stiffness
G	shear modulus for Saint-Venant torsion	torsional stiffness
A	cross-sectional area	axial stiffness E*A/L
Iy	second moment of area about local 2	bending response in local 1-3 plane, coupled w and th2 DOFs
Iz	second moment of area about local 3	bending response in local 1-2 plane, coupled v and th3 DOFs
J	Saint-Venant torsion constant	torsional stiffness G*J/L

The I/O Definition Agent owns the final property source and ordering. This formulation requires only that the implementation receives these six physical quantities with the meanings above.

Strong Form and Boundary Conditions

Let x be the coordinate along the undeformed beam centerline, 0 <= x <= L.

For constant prismatic section properties and no element-generated distributed loads, the internal generalized resultants are:

N  = E*A  * du/dx
T  = G*J  * dth1/dx
M2 = E*Iy * dth2/dx
M3 = E*Iz * dth3/dx

The first-scope homogeneous field equilibrium is:

dN/dx    = 0
dT/dx    = 0
d2M2/dx2 = 0 in the work-conjugate Euler-Bernoulli bending statement for w/th2
d2M3/dx2 = 0 in the work-conjugate Euler-Bernoulli bending statement for v/th3

Essential boundary conditions are prescribed nodal translations and rotations. Natural boundary terms are nodal axial force, shear forces, torsional moment, and bending moments. Applied nodal forces and moments are supplied by Abaqus outside the first-scope UEL.

Weak or Variational Form

The element internal virtual work is:

delta W_int =
  integral_0^L [
    delta eps  * E*A  * eps
  + delta kap1 * G*J  * kap1
  + delta kap2 * E*Iy * kap2
  + delta kap3 * E*Iz * kap3
  ] dx

with generalized strains:

eps  = du/dx
kap1 = dth1/dx
kap2 = dth2/dx = -d2w/dx2
kap3 = dth3/dx =  d2v/dx2

The first-scope element does not generate external virtual work internally:

delta W_ext_element = 0

Therefore the element internal force is derived entirely from delta W_int, and solver-applied nodal loads remain outside this formulation.

Discretization and Kinematics

Use xi = x/L, where 0 <= xi <= 1.

Linear axial and torsional interpolation:

N1 = 1 - xi
N2 = xi

u(x)   = N1*u1     + N2*u2
th1(x) = N1*th1_1  + N2*th1_2

Cubic Hermite interpolation for bending about local 3, using th3 = dv/dx:

H1 = 1 - 3*xi^2 + 2*xi^3
H2 = L*(xi - 2*xi^2 + xi^3)
H3 = 3*xi^2 - 2*xi^3
H4 = L*(-xi^2 + xi^3)

v(x) = H1*v1 + H2*th3_1 + H3*v2 + H4*th3_2

Cubic Hermite interpolation for bending about local 2, using th2 = -dw/dx:

w(x) = H1*w1 - H2*th2_1 + H3*w2 - H4*th2_2

Expected shape-function checks:

• Axial and torsion interpolation: N1 + N2 = 1.
• Hermite displacement interpolation: H1(0)=1, H3(1)=1, and the opposite endpoint displacement shape values are zero.
• Hermite slope interpolation recovers dv/dx = th3 and dw/dx = -th2 at the corresponding nodes.

Local Element Stiffness Matrix

Define:

L  = element length
EA = E*A
GJ = G*J
B2 = E*Iy
B3 = E*Iz

a  = EA/L
t  = GJ/L
b2_12 = 12*B2/L^3
b2_6  =  6*B2/L^2
b2_4  =  4*B2/L
b2_2  =  2*B2/L
b3_12 = 12*B3/L^3
b3_6  =  6*B3/L^2
b3_4  =  4*B3/L
b3_2  =  2*B3/L

In the local DOF order

[u1, v1, w1, th1_1, th2_1, th3_1, u2, v2, w2, th1_2, th2_2, th3_2]

the local stiffness matrix is:

k_local =
[
 [ a,      0,      0,   0,      0,      0, -a,      0,      0,   0,      0,      0],
 [ 0,  b3_12,      0,   0,      0,  b3_6,  0, -b3_12,      0,   0,      0,  b3_6],
 [ 0,      0,  b2_12,   0, -b2_6,      0,  0,      0, -b2_12,  0, -b2_6,      0],
 [ 0,      0,      0,   t,      0,      0,  0,      0,      0,  -t,      0,      0],
 [ 0,      0, -b2_6,   0,  b2_4,      0,  0,      0,  b2_6,   0,  b2_2,      0],
 [ 0,  b3_6,      0,   0,      0,  b3_4,  0, -b3_6,      0,   0,      0,  b3_2],
 [-a,      0,      0,   0,      0,      0,  a,      0,      0,   0,      0,      0],
 [ 0, -b3_12,     0,   0,      0, -b3_6,  0,  b3_12,      0,   0,      0, -b3_6],
 [ 0,      0, -b2_12,  0,  b2_6,      0,  0,      0,  b2_12,  0,  b2_6,      0],
 [ 0,      0,      0,  -t,      0,      0,  0,      0,      0,   t,      0,      0],
 [ 0,      0, -b2_6,   0,  b2_2,      0,  0,      0,  b2_6,   0,  b2_4,      0],
 [ 0,  b3_6,      0,   0,      0,  b3_2,  0, -b3_6,      0,   0,      0,  b3_4]
]

Subblock interpretation:

• Axial u1/u2 block: EA/L * [[1, -1], [-1, 1]].
• Torsion th1_1/th1_2 block: GJ/L * [[1, -1], [-1, 1]].
• Bending about local 2 uses E*Iy and couples w with th2.
• Bending about local 3 uses E*Iz and couples v with th3.

For a valid free element, k_local is symmetric positive semidefinite with six rigid-body modes and six deformational modes.

Local Frame and Transformation to Global Coordinates

Let the original global nodal coordinates be X1 and X2.

x21 = X2 - X1
L   = norm(x21)
e1  = x21 / L

Let a_ref be the formulation-level orientation reference vector expressed in global coordinates. The I/O Definition Agent will decide whether a_ref comes from properties, an additional point, or another approved Abaqus-facing source. The reference must not be parallel to e1.

Project the reference onto the plane normal to the beam axis:

a_perp = a_ref - dot(a_ref, e1)*e1
e2     = a_perp / norm(a_perp)
e3     = cross(e1, e2)

This gives a right-handed frame:

cross(e1, e2) = e3
det(R)        = +1

Define the 3-by-3 direction-cosine matrix R with local axes as rows:

R = [
  e1^T
  e2^T
  e3^T
]

For any global vector g, the local vector is:

l = R*g

The 12-by-12 element transformation is block diagonal:

T = blockdiag(R, R, R, R)

using the block order: node 1 translations, node 1 rotations, node 2 translations, node 2 rotations.

The local and global element vectors are related by:

q_local = T*q_global

Virtual work gives:

f_internal_local  = k_local*q_local
f_internal_global = T^T*f_internal_local
K_global          = T^T*k_local*T

AMATRX must receive K_global, not k_local.

Element Equations and Residual Convention

The first-scope element has no element-generated external load vector:

p_element = 0

Internal force:

f_int = K_global*q_global

Using the Abaqus research fact that the static residual is external minus internal, the formulation-level residual contribution is:

r_element = p_element - f_int = -K_global*q_global

Therefore the first-scope RHS mathematical contribution is:

RHS = -K_global*U

The I/O Definition Agent must confirm the exact Abaqus fill convention and no-Abaqus sign test before implementation. If the interface contract later chooses a different array-fill convention, it must preserve the same physical internal force f_int = K_global*U and document the sign mapping explicitly.

AMATRX and RHS Contribution Rules by Supported LFLAGS

This formulation supports only small-displacement static requests. Exact procedure-code ownership remains with the I/O Definition Agent, but the contribution rules are:

request condition	formulation contribution
small-displacement static and LFLAGS(3)=1	provide AMATRX = K_global and RHS = -K_global*U
small-displacement static and LFLAGS(3)=2	provide AMATRX = K_global; no residual contribution is mathematically requested
small-displacement static and LFLAGS(3)=5	provide RHS = -K_global*U; no stiffness contribution is mathematically requested
large-displacement request	out of scope; interface must reject or report unsupported behavior
mass, damping, perturbation-only, dynamics, or other non-static requests	out of scope; interface must reject or return a documented zero contribution with diagnostics

The element stiffness is symmetric. UNSYMM is not required by this formulation unless the interface contract introduces a non-symmetric behavior in a later approved scope.

Numerical Tolerances and Validity Checks

The first implementation should use double precision. These formulation tolerances are defaults for no-Abaqus tests and numerical review; the Numerical Review Agent may tighten them with scale-aware justification.

check	default rule
coordinate finiteness	reject if any coordinate used by the element is NaN or infinite
length scale	coord_scale = max(1.0, norm(X1), norm(X2))
zero or near-zero length	reject if L <= 1.0e-12*coord_scale
orientation vector finiteness	reject if any a_ref component is NaN or infinite
zero orientation vector	reject if norm(a_ref) <= 1.0e-12 when a_ref is supplied as a direction; if supplied as a point offset, apply the length scale rule
near-parallel orientation	reject if norm(a_perp)/norm(a_ref) <= 1.0e-8
local frame orthonormality	require maxabs(R*R^T - I) <= 1.0e-12 after construction
right-handed frame	require det(R) > 0 and preferably abs(det(R) - 1) <= 1.0e-12
physical properties	reject if E, G, A, Iy, Iz, or J is NaN, infinite, zero, or negative
matrix finiteness	reject or report invalid if any K_global or residual component is NaN or infinite
matrix symmetry	require maxabs(K_global - K_global^T) <= 1.0e-12*max(1.0, maxabs(K_global))

No positive lower bound is imposed on physical properties beyond > 0 because units and model scale are user-controlled. Conditioning risks from very small positive values must be caught by no-Abaqus tests and numerical review.

State Variables, Energy, and Output Recovery

State Variables

No persistent material or history state is required for the first-scope linear elastic beam. SVARS is not used by the formulation.

Optional diagnostic storage of local end forces or section quantities in SVARS is not approved by this formulation. It requires an explicit interface contract and reference extraction plan.

Energy

The elastic strain energy associated with the formulation is:

Ue = 0.5*q_global^T*K_global*q_global
   = 0.5*q_local^T*k_local*q_local

Populating Abaqus ENERGY is not required in the first scope. If a later interface contract requires energy output, it should use this expression for elastic strain energy and define the exact ENERGY slot.

Output Recovery

Quantities available from the formulation:

q_local          = T*q_global
f_local          = k_local*q_local
f_internal_global = K_global*q_global
r_element        = -f_internal_global
Ue               = 0.5*q_global^T*K_global*q_global

Local element end-force component order follows the local DOF order:

[N1, V2_1, V3_1, T1, M2_1, M3_1, N2, V2_2, V3_2, T2, M2_2, M3_2]

Initial external Abaqus CSV validation should compare:

• nodal displacements U in global coordinates at selected nodes
• support reactions RF in global coordinates at selected constrained nodes
• coordinate-system labels and declared units for every compared CSV value

External element force or section-force CSV comparison is deferred until the interface defines an approved output route, such as diagnostic SVARS.

Stress and strain recovery is not part of the first scope.

Mapping and Numerical Integration

The element is a straight two-node member, so the centerline mapping is:

X(xi) = (1 - xi)*X1 + xi*X2
dX/dx = e1
dx/dxi = L

The stiffness matrix is evaluated analytically from the closed-form axial, torsion, and Hermite bending expressions. No Gauss integration loop is required for the first implementation.

If a later implementation chooses numerical integration internally, it must reproduce the analytical matrix above within the no-Abaqus stiffness tolerances. Reduced integration is not needed and must not be introduced as an hourglass-prone substitute for the closed-form matrix.

Algorithm Pseudocode

input:
  X1, X2
  orientation reference a_ref
  E, G, A, Iy, Iz, J
  global element displacement vector U
  supported static request flag

validate:
  finite coordinates and properties
  positive E, G, A, Iy, Iz, J
  valid length L
  valid nonparallel orientation reference

construct local frame:
  e1 = (X2 - X1)/L
  e2 = normalized projection of a_ref normal to e1
  e3 = cross(e1, e2)
  R  = rows(e1, e2, e3)
  T  = blockdiag(R, R, R, R)

assemble local stiffness:
  compute EA/L, GJ/L, E*Iy/L terms, E*Iz/L terms
  fill symmetric 12 by 12 k_local in the documented local DOF order

transform:
  K_global = transpose(T) * k_local * T

compute requested contributions:
  if stiffness requested:
    AMATRX = K_global
  if residual requested:
    RHS = -K_global * U

optional recovery for tests:
  q_local = T * U
  f_local = k_local * q_local
  strain_energy = 0.5 * transpose(U) * K_global * U

Numerical Risks

risk	expected check or mitigation
rigid_body_modes	free-element K_global must have six rigid-body modes; no-Abaqus tests should check near-zero residual for rigid translations and rotations
matrix symmetry	K_global must satisfy the symmetry tolerance; asymmetry indicates transform, fill, or indexing error
positive definiteness	unconstrained element is positive semidefinite; constrained benchmark models should be positive definite after essential BCs remove rigid modes
local axis singularity	zero length and near-parallel orientation references must be rejected before normalization
axis/property swap	separate Iy and Iz bending tests are required to catch local 2/local 3 mapping mistakes
residual sign error	no-Abaqus test must verify RHS = -K_global*U under the approved interface convention
conditioning	very short elements, very long elements, or extreme property ratios can produce ill-conditioned matrices; numerical review should define scale-aware diagnostics
shear locking	not applicable to the element stiffness because Euler-Bernoulli shear deformation is not represented; applicability to deep beams remains out of scope
hourglass	not applicable because the first formulation uses a closed-form stiffness matrix rather than reduced integration
volumetric locking	not applicable
section offset and warping	out of scope; tests must not expect these effects

Verification-Relevant Invariants

No-Abaqus verification should be able to check:

• identity-orientation K_global = k_local
• rotated-orientation K_global = T^T*k_local*T
• maxabs(K_global - K_global^T) within tolerance
• six free-element rigid-body modes
• axial response proportional to E*A/L
• torsional response proportional to G*J/L
• local 2 bending response proportional to E*Iy
• local 3 bending response proportional to E*Iz
• residual consistency with RHS = -K_global*U

External Abaqus comparison should initially check nodal displacement and support reaction CSVs only, with explicit units and coordinate-system labels.

Open Issues and Downstream Handoff

Numerical Review Agent

• Review the local stiffness signs for th2 = -dw/dx and th3 = dv/dx.
• Review rigid-body mode checks, matrix rank expectations, and symmetry tolerance.
• Review whether 1.0e-8 near-parallel orientation rejection is appropriate for double-precision production use.
• Review conditioning risks for extreme L, Iy/Iz, J, and property-ratio values.

I/O Definition Agent

• Define the exact Abaqus-facing source and order for E, G, A, Iy, Iz, J, and a_ref.
• Define whether invalid inputs are reported by PNEWDT, message output, hard stop, or another approved Abaqus-compatible mechanism.
• Confirm the exact RHS fill convention and static LFLAGS behavior before implementation.
• Decide whether ENERGY remains unused or receives elastic strain energy in a later approved scope.
• Keep SVARS unused unless a diagnostic output contract is explicitly approved.

Reference Model Agent

• Include no-Abaqus tests for axial, torsion, local 2 bending, local 3 bending, transformation, rigid-body modes, matrix symmetry, invalid orientation, invalid length, invalid properties, and residual sign.
• Keep external Abaqus evidence as user-generated model.inp, ODB-extracted nodal U/RF CSV files, log tails, and metadata under references/uel-3d-euler-beam/<model-id>/.
• Do not require element force CSV comparison until an output recovery path is approved.

Implementation Planning Agent

• Implement only after no-Abaqus RED evidence exists.
• Keep the Abaqus UEL wrapper thin and isolate the beam matrix, transform, validation, and residual calculations in testable routines if the approved interface permits it.
• Preserve the local DOF order and transformation convention exactly as documented here unless Numerical Review and I/O Definition update this formulation.