FEM Formulation

This section describes the finite element formulation used by the Armatura solver. The implementation follows the direct stiffness method as presented in Kassimali (2015) and the CALFEM 3.6 toolbox from Lund University.

Readers who want to build their own understanding of the underlying theory should consult the primary references directly — this page provides orientation, not a substitute for the textbooks.

Primary references

Shortcode	Source	Used for
MAS	Kassimali, A. (2015). Matrix Analysis of Structures, 3rd ed. Cengage Learning.	Stiffness method, DOF numbering, transformations, element derivations, distributed load integration
CALFEM	Austrell et al. (2022). CALFEM — A Finite Element Toolbox, Version 3.6. Lund University.	Element formulations (beam3e), shape functions, orientation vector approach
EC3	EN 1993-1-1:2005. Design of Steel Structures — Part 1-1. CEN.	Material properties, design checks
SN003a	SN003a-EN-EU. Lateral Torsional Buckling. SCI.	Elastic critical moment Mcr

Direct stiffness method

The solver uses the classical direct stiffness method. The steps are:

1. Compute the local stiffness matrix for each element
2. Transform each local stiffness matrix to global coordinates
3. Assemble the global stiffness matrix by scattering element contributions
4. Apply boundary conditions (support constraints)
5. Solve the system Ku = F for unknown displacements
6. Back-substitute to recover reactions and member forces

This procedure is described in [MAS] Chapter 2 (§2.1–2.6). The implementation follows Kassimali's DOF numbering and sign conventions throughout.

Euler–Bernoulli beam element

All beam members use the 3D Euler–Bernoulli (EB) beam element with 6 degrees of freedom per node (3 translations + 3 rotations), giving a 12×12 element stiffness matrix.

The EB formulation assumes:

• Plane sections remain plane and normal to the neutral axis after deformation
• Shear deformation is neglected (valid for slender members where L/h > 10)
• Small displacements and rotations (linear elastic analysis)

The local stiffness matrix is assembled from four decoupled sub-problems: axial (EA/L), torsion (GJ/L), bending in the x-y plane (EIz), and bending in the x-z plane (EIy). The full 12×12 matrix and its derivation are given in [MAS] §8.3, Eq. 8.39 (p. 469) and [CALFEM] §5.6.

The bending stiffness terms use cubic Hermite interpolation functions. These shape functions are exact for prismatic elements under end loads — meaning the FEM solution at the nodes matches the analytical solution exactly for nodal loading. For loads between nodes (distributed loads), equivalent nodal forces are used.

When EB is appropriate

The EB formulation is accurate for typical steel beams and columns where L/h > 10. This covers the vast majority of structural steel members in practice. For very short, deep members (L/h < 5), Timoshenko beam theory would give better results by including shear deformation — this is not currently implemented.

Transformation matrix

Each element's local stiffness matrix must be transformed to global coordinates before assembly. The transformation uses a 12×12 block-diagonal matrix built from a 3×3 direction cosine matrix (λ).

The implementation uses a CALFEM-style orientation vector (eo) to define the element's local coordinate system, rather than a rotation angle γ. This approach uses cross products to derive the local axes:

1. Local x — unit vector from start node to end node
2. Local z — component of the orientation vector perpendicular to local x (Gram–Schmidt)
3. Local y — cross product of local z and local x

This is the same algorithm used in CALFEM's beam3e function ([CALFEM] §5.6). It handles all member orientations including vertical members without special cases — the orientation vector simply must not be parallel to the member axis.

For non-vertical members, the default orientation vector is global Z = (0, 0, 1). For vertical members (where global Z is parallel to the member axis), the default falls back to global X = (1, 0, 0). Users can override this per member via Member.OrientationVector.

The transformation is described in [MAS] §8.3 (Eq. 8.35–8.38). The orientation vector approach is documented in [CALFEM] §5.6 (the eo parameter of beam3e).

A note on earlier approaches

An earlier version of the solver used a gamma-angle formulation with explicit special cases for vertical members. This caused subtle bugs where section Ixx/Iyy values were mapped to the wrong bending plane for non-horizontal members. The orientation vector approach eliminated these bugs by making the axis mapping explicit and consistent for all orientations. See the Sign Conventions page for details on how section properties map to local axes.

Section property mapping

The element stiffness matrix requires Iy (bending about local y) and Iz (bending about local z) in element-local coordinates.

The mapping from section properties is fixed:

• Section.Ixx → Iy (resists bending in the x-z plane, causing z-deflection)
• Section.Iyy → Iz (resists bending in the x-y plane, causing y-deflection)

This is the same convention as CALFEM: the user provides Iy and Iz already resolved to local axes. In Armatura, the orientation vector ensures that the local axes align with the section's physical axes, so the mapping is always Ixx → Iy, Iyy → Iz — no conditional swapping.

See [MAS] §8.3, p. 469 (note after Eq. 8.39) for the relationship between section inertias and element bending stiffness.

Assembly and boundary conditions

The global stiffness matrix is assembled using the standard scatter operation: each element's transformed stiffness contribution is added to the global matrix at the DOF positions corresponding to its nodes. The DOF numbering scheme follows [MAS] §8.3 (p. 462): six DOFs per node in the order [Ux, Uy, Uz, Rx, Ry, Rz].

Boundary conditions are applied by modifying the corresponding rows and columns of the global stiffness matrix. Fixed DOFs are constrained by zeroing the off-diagonal terms and setting the diagonal to 1 (penalty method), with the corresponding force vector entries set to the prescribed displacement (zero for standard supports).

Equivalent nodal forces

Distributed loads acting along members are converted to equivalent nodal forces using Gauss quadrature over the cubic Hermite shape functions. This is the consistent load approach from finite element theory.

Gauss quadrature is the production approach because it handles any load distribution (uniform, linearly varying, or arbitrary) without member subdivision. The number of Gauss points is chosen to integrate the load–shape-function product exactly for the polynomial order involved: 2 points suffice for uniform loads on prismatic elements.

The equivalent nodal force formulas for standard load types (uniform, triangular, trapezoidal) are given in [MAS] Chapter 6 and on the inside front cover table.

Point loads are converted to equivalent nodal forces using closed-form fixed-end force expressions from [MAS] §5.3 (inside front cover table). For a transverse force P at distance a from the start of an element of length L (with b = L - a): V1 = Pb²(3a+b)/L³, M1 = Pab²/L², V2 = Pa²(a+3b)/L³, M2 = -Pa²b/L². Axial and torsional components use linear interpolation. Internal force recovery adds a piecewise span correction at the load point: a step in shear and a slope change in moment ([MAS] §6.4).

Internal force and displacement recovery

After solving for nodal displacements, internal forces and displacements at intermediate points along a member are recovered using shape function interpolation plus a span load correction (particular solution). The correction accounts for the difference between the cubic shape function approximation and the true higher-order displacement field under distributed loading.

For uniform loads, the closed-form particular solution from [MAS] §6.4 is used directly. This avoids subdivision while giving exact results at any point along the member.

Results are reported at configurable stations along each member (default: 11 stations from 0.0 to 1.0 in increments of 0.1).

Solver

The global system Ku = F is solved using LU factorization via MathNet.Numerics with LAPACK native binaries. For the problem sizes typical of frame analysis (up to a few thousand DOFs), this direct solver approach is efficient and numerically stable.

For larger problems (10,000+ DOFs with sparse matrices), an iterative solver or sparse direct solver would be beneficial. This is not currently implemented but would be a natural extension — the assembly and solution phases are already separated.

Limitations of the current formulation

The solver currently implements first-order linear elastic analysis only. This means:

• No second-order effects (P-Δ or P-δ)
• No geometric nonlinearity
• No material nonlinearity (plasticity)
• No Timoshenko shear deformation

Second-order effects are important for slender steel frames under compression. The αcr check (EN 1993-1-1 §5.2.1) is implemented as a screening tool to warn when first-order analysis may be unconservative. Linear eigenvalue buckling analysis and modal analysis are planned.

Rigid links and attached bodies

Attached bodies are connected to the structural frame by rigid links and a rigid spine chain. The rigid links are implemented through a penalty-stiffness material model:

• E_rigid = 1e8 * E_reference

This follows the standard penalty-method approach for near-constraint enforcement (Cook et al., 2002, §9.5). Loads from attached bodies are applied at the COG node. The solver then transfers and distributes those loads through the rigid chain, capturing eccentricity without manual tributary splitting.

Where attachment points fall inside members, analysis-generated internal nodes can be inserted at extraction/preprocessing time. These nodes are analysis artifacts and do not need to be user-authored editor nodes.

Further reading

For readers wanting to understand the theory in depth:

• [MAS] Chapters 2, 3, 6, and 8 cover the complete direct stiffness method for 3D frames
• [CALFEM] Sections 5.4–5.6 cover the CALFEM beam element implementations
• Cook, Malkus, Plesha & Witt (2002). Concepts and Applications of Finite Element Analysis, 4th ed. — broader FEM context beyond frames