FESADev/docs/formulations/euler-beam-3d-formulation.md

3D Euler-Bernoulli Beam Kernel Formulation

Metadata

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

Scope and Assumptions

• analysis_type: linear static kernel only.
• element_type: two-node straight prismatic 3D Euler-Bernoulli beam.
• deformation: small displacement and small rotation.
• linearity: linear elastic.
• material_model_boundary: constant E and G.
• section_model_boundary: constant A, J, Iy, Iz.
• coordinate_system: local Cartesian beam basis and global Cartesian basis.
• units: user-consistent.

Primary Variables and DOFs

• nodal_variables: displacement vector and rotation vector at each node.
• dof_ordering: [u1, v1, w1, rx1, ry1, rz1, u2, v2, w2, rx2, ry2, rz2].
• local_axis_convention:
  • local x is the normalized vector from node 1 to node 2.
  • local y is the normalized projection of the user orientation vector onto the plane normal to local x.
  • local z = x cross y.
• sign_convention: internal local end force vector is f_local = K_local * u_local in the same DOF order.

Strong Form and Boundary Conditions

The kernel is an element-level direct-stiffness implementation of the linear Euler-Bernoulli frame equations. For a prismatic member in local coordinate x:

• axial equilibrium uses constant axial stiffness EA.
• torsional equilibrium uses constant torsional stiffness GJ.
• bending in local x-y uses flexural rigidity EIz.
• bending in local x-z uses flexural rigidity EIy.

Essential and natural boundary conditions are not applied by this kernel. They are downstream solver assembly concerns.

Weak or Variational Form

The element stiffness represents the second variation of strain energy:

U = 1/2 integral_0^L [
      EA * (du/dx)^2
    + GJ * (drx/dx)^2
    + EIz * (d2v/dx2)^2
    + EIy * (d2w/dx2)^2
] dx

The element internal force is the derivative of this energy with respect to the nodal DOFs:

f_local = K_local * u_local

External load vectors are out of scope for this increment.

Discretization

• axial interpolation: linear two-node interpolation for u.
• torsion interpolation: linear two-node interpolation for rx.
• bending interpolation: cubic Hermite interpolation for each transverse displacement and corresponding nodal rotation.
• nodal_layout: node 1 at local coordinate x=0, node 2 at x=L.
• partition_of_unity_check: axial and torsion linear shape functions sum to one.
• kronecker_delta_check: axial, torsion, and Hermite bending shape functions recover the nodal displacement/rotation DOFs.

Kinematics

• axial strain: epsilon_x = du/dx.
• torsion rate: kappa_x = drx/dx.
• bending curvature for local x-y bending: curvature associated with v and rotations rz.
• bending curvature for local x-z bending: curvature associated with w and rotations ry.
• strain_measure: infinitesimal strain and small rotation curvature measures only.

Constitutive Contract

• axial force result: N = EA * du/dx.
• torsional moment result: Mx = GJ * drx/dx.
• bending moment about local z: uses EIz.
• bending moment about local y: uses EIy.
• material_state_variables: none.

Element Equations

Let:

a  = E*A/L
t  = G*J/L
by = E*Iy
bz = E*Iz
cy1 = 12*by/L^3
cy2 =  6*by/L^2
cy3 =  4*by/L
cy4 =  2*by/L
cz1 = 12*bz/L^3
cz2 =  6*bz/L^2
cz3 =  4*bz/L
cz4 =  2*bz/L

K_local is a 12x12 symmetric matrix initialized to zero, with these nonzero entries before symmetric mirroring:

K(0,0)= a       K(0,6)=-a       K(6,6)= a

K(3,3)= t       K(3,9)=-t       K(9,9)= t

K(1,1)= cz1     K(1,5)= cz2     K(1,7)=-cz1    K(1,11)= cz2
K(5,5)= cz3     K(5,7)=-cz2     K(5,11)= cz4
K(7,7)= cz1     K(7,11)=-cz2    K(11,11)= cz3

K(2,2)= cy1     K(2,4)=-cy2     K(2,8)=-cy1    K(2,10)=-cy2
K(4,4)= cy3     K(4,8)= cy2     K(4,10)= cy4
K(8,8)= cy1     K(8,10)= cy2    K(10,10)= cy3

The implementation must mirror upper-triangular entries to the lower triangle.

Mapping and Numerical Integration

• reference_coordinates: local beam coordinate x in [0,L].
• isoparametric_mapping: straight member mapping only.
• jacobian: dx/dxi = L/2 if a parent coordinate is introduced; the direct closed-form stiffness does not require runtime quadrature.
• determinant_checks: L > 0 and finite.
• gauss_points_and_weights: not used by the closed-form kernel.
• integration_policy: analytical closed-form.

Transformation

Let R be the 3x3 matrix whose rows are local basis vectors expressed in global coordinates:

R = [ x_local^T
      y_local^T
      z_local^T ]

For each node, the same R maps global translational and rotational components to local components. The 12x12 transform T is block diagonal:

u_local = T * u_global
K_global = T^T * K_local * T

Global end-force recovery uses the same convention:

u_local = T * u_global
f_local = K_local * u_local
f_global = T^T * f_local

Output Recovery

• displacement: input only for this kernel.
• reaction: not recovered by this kernel; downstream assembly/constraints own reactions.
• element_force: local and global element end forces recover from stiffness times displacement.
• strain: out of scope.
• stress: out of scope.
• nodal_extrapolation: not applicable.

Invalid Input Handling

• length <= 0, nonfinite length, or near-zero geometry length must throw std::invalid_argument.
• nonpositive or nonfinite E, G, A, J, Iy, or Iz must throw std::invalid_argument.
• zero orientation vector or orientation parallel to local x must throw std::invalid_argument.

Unit Test Tolerance

• representative coefficient checks: absolute tolerance 1.0e-10 for the planned deterministic values.
• symmetry checks: absolute tolerance 1.0e-10.
• rigid mode force norm checks: absolute tolerance 1.0e-9 for simple test magnitudes.

Algorithm Pseudocode

validate length and section constants
compute named stiffness coefficients
initialize 12x12 local matrix to zero
fill upper triangular axial, torsion, and bending terms
mirror upper triangular terms

for global operations:
  compute normalized local x from node2 - node1
  project orientation onto plane normal to local x and normalize as local y
  compute local z = x cross y
  build block diagonal T from R
  compute K_global = T^T * K_local * T
  recover forces with u_local = T*u_global, f_local = K_local*u_local, f_global = T^T*f_local

Numerical Risks

• rigid_body_modes: the unconstrained local/global stiffness should produce near-zero force for six rigid body modes; tests should include at least rigid translation in this increment.
• patch_test: future solver-level patch tests require parser/assembly integration.
• symmetry: local and global matrices must remain symmetric.
• positive_definiteness: constrained systems may become positive definite; the isolated element is positive semidefinite.
• hourglass: not applicable to this closed-form beam kernel.
• shear_locking: excluded by Euler-Bernoulli assumptions and no shear deformation terms.
• volumetric_locking: not applicable.
• distortion: curved or offset geometry is out of scope.
• singular_jacobian: represented by zero or near-zero length.

Open Issues and Downstream Handoff

Numerical Review Agent

• Confirm signs for the local w and ry bending block.
• Confirm the transformation convention and invalid orientation handling are implementation-ready.

I/O Definition Agent

• Define future orientation input mapping and unsupported cases without claiming parser support in this kernel increment.

Reference Model Agent

• Define future axial, torsion, bending, and skew-oriented reference models without generating artifacts.

Implementation Planning Agent

• Transcribe the matrix entries and transformation convention into C++ tests before production implementation.