# 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.

```text
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:

```text
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:

```text
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:

```text
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:

```text
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:

```text
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:

```text
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:

```text
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`:

```text
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`:

```text
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:

```text
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

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

the local stiffness matrix is:

```text
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`.

```text
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:

```text
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:

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

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

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

For any global vector `g`, the local vector is:

```text
l = R*g
```

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

```text
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:

```text
q_local = T*q_global
```

Virtual work gives:

```text
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:

```text
p_element = 0
```

Internal force:

```text
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:

```text
r_element = p_element - f_int = -K_global*q_global
```

Therefore the first-scope `RHS` mathematical contribution is:

```text
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:

```text
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:

```text
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:

```text
[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:

```text
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

```text
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.