Finite Element Analysis (FEA)

The conceptual foundations of the finite element method (FEM) trace back to the early 1940s. In 1941, Alexander Hrennikoff, a Russian-Canadian structural engineer at the University of British Columbia, published a paper describing the approximation of continuous membrane and plate structures using a lattice framework of discrete elements. In 1943, the mathematician Richard Courant introduced the use of triangular and rectangular sub-domains (effectively early 'finite elements') for solving the torsion problem of a prismatic bar, applying variational calculus to approximate solutions of partial differential equations. However, it was the engineering community that developed FEM into a practical computational tool. In 1956, M.J. Turner, R.W. Clough, H.C. Martin, and L.J. Topp at Boeing published a seminal paper presenting the direct stiffness method for structural analysis using discrete triangular elements. Ray W. Clough formally coined the term 'finite element method' in a 1960 conference paper titled 'The Finite Element Method in Plane Stress Analysis,' presented at the 2nd ASCE Conference on Electronic Computation. Through the 1960s and 1970s, Olgierd Zienkiewicz at Swansea University, along with John H. Argyris at Stuttgart, developed the mathematical and computational frameworks that made FEM widely applicable across civil, mechanical, and aerospace engineering. By the 1980s and 1990s, advances in computing power and the development of commercial FEA software packages such as ANSYS, Abaqus, and Strand7 made the method accessible to a much broader research community, including biologists and paleontologists.

The transfer of FEA from engineering to biological research began in the orthopedic and dental biomechanics fields, where FE models were used to study the mechanical behavior of bones, implants, and teeth under functional loading. By the late 1990s, the method began to be applied to questions in vertebrate functional morphology. The landmark application of FEA to paleontology came in 2001, when Emily J. Rayfield and colleagues published a study in Nature analyzing the skull of the large Jurassic theropod dinosaur Allosaurus fragilis. This study generated one of the most geometrically complex FE models of any organism at that time, demonstrating that the Allosaurus skull could withstand forces substantially greater than its estimated muscle-driven bite force — suggesting the skull was structurally 'overbuilt.' In 2004, Rayfield followed this with an analysis of the Tyrannosaurus rex skull, testing the 'puncture-pull' feeding hypothesis by simulating biting and tearing loads. The results showed the T. rex cranium was equally adapted to resist both types of loading, and that compressive and shear stresses concentrated in the nasals rather than the fronto-parietal region (as in Allosaurus), providing a biomechanical rationale for the robusticity commonly observed in tyrannosaurid nasals. The maxilla-jugal suture was found to provide a tensile shock-absorbing function, reducing localized tension but slightly weakening the skull overall. Rayfield's 2007 review article in the Annual Review of Earth and Planetary Sciences provided a comprehensive framework for FEA in paleontology and became a foundational reference for subsequent studies.

A typical paleontological FEA workflow involves several stages. First, the three-dimensional geometry of a fossil structure is acquired, usually through CT scanning (micro-CT or medical CT), though photogrammetry and surface laser scanning are also used. The resulting volumetric data are segmented using software such as Avizo, Mimics, or open-source alternatives like 3D Slicer, isolating the bony structure of interest from surrounding matrix. The segmented surface mesh undergoes cleaning and repair to remove artifacts such as holes, non-manifold edges, overlapping faces, and excessively heterogeneous element sizes. The cleaned surface mesh is then converted into a volumetric mesh typically composed of four-noded tetrahedral elements (tet4), though higher-order elements and other geometries (shell, beam, plate) can be used for specific applications.

Material properties must be assigned to the mesh elements. For most paleontological studies, bone is modeled as a homogeneous, isotropic, linearly elastic material, with typical values of Young's modulus (E) around 17,000–20,000 MPa and Poisson's ratio (ν) around 0.3, derived from experimental testing of extant archosaur bone. Some studies assign heterogeneous or anisotropic properties, and recent work has explored the effects of modeling tooth enamel (E ≈ 50,000–80,000 MPa) separately from dentine (E ≈ 18,000 MPa).

Boundary conditions are then defined: muscle forces are applied to attachment areas on the model surface (estimated from dry skull methods, muscle reconstruction, or physiological cross-sectional area calculations), and constraints (fixed displacement points) are placed at biologically relevant locations such as the temporomandibular joint or tooth–food contact points. Loading scenarios simulate different behaviors such as unilateral or bilateral biting at various tooth positions. The model is then solved using a linear static solver in commercial software (ANSYS, Abaqus, Strand7) or increasingly in open-source tools (Fossils, FEBio). The output includes distributions of von Mises stress, principal stresses and strains, deformation, strain energy, and reaction forces (bite forces) at constrained nodes.

Several mechanical metrics are commonly extracted from paleontological FE models. Von Mises stress is the most widely reported, as it provides a scalar measure combining the effects of all stress components, useful for predicting yielding in ductile materials like bone. Principal stresses (maximum and minimum) reveal the orientation and magnitude of tensile versus compressive loading. Strain energy, the total work done in deforming the structure, has been used as a measure of structural efficiency — a skull that converts more input muscle force into output bite force with less strain energy is considered more mechanically efficient. Mechanical efficiency, defined as the ratio of bite reaction force to total input muscle force, allows direct comparison of feeding performance between taxa.

A critical challenge in paleontological FEA is model validation, since it is impossible to experimentally test the mechanical behavior of extinct organisms. Validation studies have been conducted using extant analogues. Strait et al. (2005) compared FE predictions against in vivo strain gauge data from Macaca (rhesus macaque) skulls during biting, finding that while absolute strain magnitudes were sometimes inaccurate, general patterns of strain distribution and orientation were well predicted by the model when muscle forces and material properties were accurately represented. Similar validation work has been performed on Alligator mississippiensis mandibles (Porro et al., 2011), pig skulls (Bright and Gröning, 2011), and lizard crania (Dutel et al., 2021). These studies collectively demonstrate that while absolute magnitudes of stress and strain must be interpreted cautiously, relative comparisons of stress distribution patterns between models built using identical protocols are robust and scientifically informative.

Paleontological FE models are built upon numerous simplifications and assumptions, and a substantial body of sensitivity literature has examined their effects. Key parameters that influence model output include mesh density (coarser meshes can produce inaccurate stress estimates), material property assignments (inclusion of enamel, dentine, and sutures can affect local stress patterns and bite force estimates), muscle force magnitude and direction (including the effect of muscle wrapping versus straight-line force vectors), boundary condition placement (number and location of constrained nodes), and the inclusion or exclusion of anatomical features such as cranial sutures, the chondrocranium, periodontal ligaments, and internal anatomy. Herbst et al. (2021) demonstrated that simplifying teeth as all-bone versus modeling realistic enamel layers did not visibly change overall stress distribution patterns, but did affect absolute bite force values by up to approximately 28%. Rayfield (2019) noted that inappropriate modeling or omission of cranial sutures can lead to inaccurate stress, strain, and deformation results.

The vast majority of paleontological FEA studies use linear static analysis, which assumes small deformations, linearly elastic materials (Hooke's Law), and static equilibrium. This approach is computationally inexpensive and sufficient for many biomechanical questions about skeletal elements. However, as reviewed by Marcé-Nogué (2022), nonlinear FEA offers important capabilities that are underutilized in paleontology. Material nonlinearities allow modeling of hyperelastic soft tissues (ligaments, cartilage, periodontal ligaments) and plastic deformation of bone. Geometric nonlinearities (large displacements) enable analysis of buckling in slender structures such as long limb bones of large dinosaurs or pterosaur wing bones. Contact nonlinearities permit the simulation of articulated multi-bone systems where separation and sliding between bone surfaces is possible, such as carpal bones, auditory ossicles, or vertebral columns. These nonlinear approaches require iterative solvers (e.g., Newton-Raphson method) and substantially greater computational resources, but can address questions that linear models cannot, including retrodeformation of taphonomically distorted fossils and the structural failure of skeletal elements under extreme loading.

FEA has become an essential tool for comparative functional morphology, allowing researchers to test hypotheses about the evolutionary optimization of skeletal form. By building FE models for multiple taxa within a clade and subjecting them to standardized loading conditions, researchers can quantify biomechanical performance differences and correlate them with ecological variables such as diet, feeding strategy, or locomotor mode. Studies have examined cranial biomechanics across theropod dinosaurs (including the macroevolutionary trends leading to extreme osteophagy in tyrannosaurids), carnivoran mammals, crocodylomorphs, temnospondyl amphibians, and many other groups. Between 2005 and 2020, more than 750 publications applied FEA to biological and evolutionary questions (Tseng, 2021). However, Tseng (2021) demonstrated an important methodological caution: small taxonomic sample sizes (fewer than approximately 11–12 taxa in 3D FE studies) can produce artificially elevated correlations between biomechanical performance and feeding ecology, generating false positive results. This finding has led to recommendations for careful study design, including attention to representative taxonomic sampling and consideration of phylogenetic structure, or alternatively adopting theoretical morphospace approaches.

Paleontological FEA has traditionally relied on commercial engineering software such as Strand7, Abaqus, ANSYS, and Hypermesh. More recently, open-source alternatives have been developed specifically for biological applications. The Fossils software (Chatar et al., 2023) provides a free, effective tool for simulating muscle-driven biomechanical loading of bone, with performance comparable to or exceeding that of commercial packages. The BFEX add-on (Díaz de León Muñoz et al., 2025) integrates Fossils with the open-source 3D software Blender, streamlining the model preparation workflow. FEBio (Maas et al., 2012) offers specialized capabilities for biological finite element analysis. These developments are democratizing access to FEA, enabling researchers without expensive software licenses to conduct rigorous biomechanical analyses.

Recent methodological advances include the integration of deep learning for automated CT segmentation (Zhang et al., 2025), combination of FEA with multibody dynamics analysis (MDA) to simulate dynamic feeding sequences rather than static snapshots, and the coupling of FEA results with geometric morphometrics and phylogenetic comparative methods to study form-function evolution across entire clades. The field continues to grow rapidly, with ongoing efforts to improve model realism through nonlinear material models, multi-body contact simulations, and better characterization of soft tissue properties via phylogenetic bracketing of extant relatives.